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accurate  results  as  a  guide,  use  FORTRAN  to  obtain  speed  and  acceptable 
accuracy.  We  believe  that  this  dual  thrust  will  finally  crack  the 
"intractable"  problem  of  STO  multicenter  integrals.  Applications  will  soon 
be  made  to  real  molecules. 


PAPERS  PUBLISHED  AND  WRITTEN 


1.  H.W.  Jones,  "Lowdin  Alpha  Function,  Overlap  Integrals,  and 
Computer  Algebra",  Int.  J.  Quantum  Chem.  (1992). 

2.  H.W.  Jones,  "Benchmark  Values  for  Two-Center  Coulomb  Integrals 
over  Slater-Type  Orbitals". 

3.  H.W.  Jones,  "Analytic  Lowdin  Alpha  Function  Method  for  Two-Center 
Electron  Repulsion  Integrals  Over  Slater-Tupe  Orbitals",  J. 
Computational  Chemistry  (December,  1991). 

4.  H.W.  Jones,  "Semi-Analytical  Method  for  Four-Center  Molecular 
Integrals  over  Slater-Type  Orbitals",  Int.  J.  Quantum  Chem.  (1992). 

5.  H.W.  Jones  and  B.  Eteamdi,  "Multicenter  Molecular  Integrals  Using 
Harmonic  Expansions  of  Slater-Type  Orbitals  and  Numerical 
Integrations",  Int.  J.  Quantum  Chem.  Symp.  24,  405  (1990). 

6.  H.W.  Jones,  "Analjrtical  Method  for  Two-  and  Three-Center  Molecular 
Integrals  over  Slater-Type  Orbitals  using  Expansions  in  Spherical 
Harmonics". 


Lowdin  Alpha-Function,  Overlap  Integral,  and  Computer  Algebra 

Herbert  W.  Jones 
Physics  Department 
Florida  A&M  University 
Tallahassee,  Florida,  32307,  U.S.A. 

Abstract 

A  commerical  computer  algebra  programme,  Mathematica,  is  used  to 
generate  the  C  matrix  that  characterizes  our  implementation  of  the  Lowdin 
alpha-function  method  as  applied  to  Slater-type  orbitals.  An  example  of  a 
two-center  overlap  integral  is  done  to  show  how  the  ai-bitrary  precision  capability 
of  Mathematica  can  overcome  severe  cancellation  errors  encountered  with 
programming  in  FORTRAN.  This  strategy  is  capable  of  being  generalized  to 
other  multicenter  molecular  integrals.  Mathematica  programmes  are  included. 
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Introduction 


Efforts  continue  to  be  made  to  make  Slater-type  orbitals  (STOs)  a  viable 
alternative  to  Gaussian-type  orbitals  (GTOs)  when  dealing  with  problems  of  ab 
initio  chemistry  [1].  The  more  physical  nature  of  STOs  should  give  advantages 
over  GTOs  when  working  in  structural  chemistry  and  in  studies  of  reaction 
dynamics,  and  density  functional  theory  [2].  Other  investigators  [3]  have  recently 
made  considerable  progress  in  this  essentially  multicenter  molecular  integral 
problem. 

In  this  paper,  we  will  initiate  a  new  line  of  action  that  takes  advantage  of  the 
availability  of  commerical  computer  algebra  programmes  of  evolving  power  and 
flexibility.  We  will  repeat  our  work  on  the  L'o'wdin  alpha-function  and,  by  an 
example,  show  how  overlap  integrals  may  be  computed  to  arbitrary  accuracy. 

The  Alpha-Function  and  C,  E,  and  F  Matrices 

The  Lowdin  a-function  method  expands  displaced  orbitals  in  an  infinite 
series  of  spherical  harmonics  about  a  single-center  (the  origin)  before  the 
required  integration  is  attempted  [4].  The  functional  coefficients  of  the  spherical 

harmonics  are  designated  as  a-functions.  We  take  a  displaced  STO  to  be  centered 

at  (0,0, a)  in  its  local  cooidinate  system  (R,  0,  (J)),  and  write  it  in  terms  of  the 

original  coordinate  system  (r,0,  (j))  [5,6].  Thus, 

X  =  ARN-1  e-CR  YM  (0,  0) , 

Li 


X=  A 

*  (2L -hi)  (L  +  M)!  ■ 

1/2  « 

E 

4ji(^-hM)!  ■ 

^N-l 

4x(L-M)! 

f=M 

(2£+l)(^-M)! 

-,1/2 


(.1)M 


X  YM(0,(})), 
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where 


N+L+i  N+.t 

aNLM  Cr)  =  (2^  +1)  (I  +  M)!  I  Z  C^LM  (ij) 

2U+M)!  i=0  j=0 


X  Hij  (Ca)i-L-^-l  (^r)K-l 

and 


e'C®  [(-l)i  eC*"  -  e'C*']  ,  r<a 

Hij  = 

[(-1)'  eC“  -  e‘Co]  ,  r>a. 


A  =  (2C)N+1/2  [(2N)!]-i/2  is  the  normalization  factor,  N,  L,  and  M  are  the  quantum 
numbers  of  the  orbital,  and  C  is  the  screening  constant  or  orbital  exponent,  r  - 

^j'JK  S  '2  ,_S.  C4  A  I-*  c  S"(“ 

U;  '  '  f  / 

Most  importantly  for  our  developments  the  elements  of  the  C  matrix  are 
integers.  Originally,  it  was  obtained  by  programming  the  following  expression, 
using  FORTRAN  and  a  simple  in-house  version  of  computer  algebra  [7]: 


[(L+M)/2] 

L+M-2p 

LrfM-2p-q 

z  z 

CNLM  (ij)ai 

ri  =  Z 

Z 

Z 

i=0  j=0 

p=0 

q=0 

v=0 

[(£  -M)/2] 

£-M-2p' 

t  -M-2p'-q'  t 

t-k 

X  Z 

Z 

Z  Z 

Z 

p’=0 

q'=0 

v'=0  k=0 

k'=0 

aX  ry  (-DV-Hq’-f-p-i-p’+L  (2L-2p)!  (2i-2p*)! 

4L+  ^+p-p'  (L-p)!  p!  pM  q!  q-j  v!  v'!  (L-i-M-2p-q-v)! 

(N-L-f2p+2q+2q')! 

(^-p')!  ( ■t-M-2p'-q'-v')!  k'!  (N-L+2p+2q+2q’-k-k’)! 
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where 


X  =  N  +  L  +  2l-2p'-  2v'-  2v  -  k  -  k’ , 
y  =  2p'  +  2v  +  2v'  +  k'  , 

and 


t  =  N-L+  2p  +2q  +  2q' 


(  .Cl 


— ^ 
— ty 


In  Table  I  we  use  the  programming  language  Mathematica  [8]  to  generate  the 
polynomial  in  a  and  r  for  the  zeroth  harmonic  (£  =  0)  for  the  2p  orbital  (N=2,  L=l, 
M=0).  (For  our  example,  we  multiply  by  (-1)^  so  that  the  positive  lobe  will  face 
toward  the  origin  [9]).  This  being  accomplished,  the  next  line  selects  the 
coefficients  of  the  polynomial  to  obtain  the  C  matrix,  whose  elements  are  given  by 
C  [[  i+1,  j+1]]  which  corresponds  to  our  form  (i,j).  Our  programming 

notation  conforms  to  Mathematica  protocol  and  is  chosen  for  ease  in  typing. 
Thus,  N=nn,  L=hh,  M=mm,  ^=h,  p'=pp,  q'=qp,  v'=vp,  k'=kp.  Floor  [  ]  means 
reduce  to  integer.  The  curly  brackets  represent  the  eight  summations  and  their 
limits.  One  of  the  great  advantages  oi Mathematica  is  that  its  high  level  language 
is  remarkably  analogous  to  standard  mathematical  notation. 

Now  that  we  have  obtained  the  C  matrix,  we  will  next  generate  the  E  matrix 

[10],  which  is  a  matrix  of  coefficients  of  the  expansion  of  the  a-function  in  a  Taylor 

series  about  r=0,  that  is,  for  the  case  r<^'  In  Table  II  we  use  two  lines  to  type  in 

I  been 

the  definitions  of  the  a-function,  alpha.  The  Mathematica  command  Series  [  ] 

expands  alpha  in  a  Taylor  series  about  r=0  to  r^.  We  again  collect  these 


coefficients  and  print  out  the  E  matrix,  e  =  E^lO  (jj).  The  F  matrix  may  be 


obtained  by  expanding  the  a-function  for  r>a  in  terms  of  a,  in  a  similar  manner, 


Overlap  Integi'al 

We  will  find  the  value  of  an  overlap  integral  given  by  Bhattacharrya  and 
Dhabal  [11].  However,  for  convenience,  our  value  will  be  positive  to  conform  to  the 
convention  of  Mulliken,  et  al.  [9].  The  definition  of  the  overlap  integral  is 

S  =  /  where 

=Aa  r^’-l  e‘C  r  Yj^,  (0,  <t>),  and  for  the  displaced  orbital  X|j=  Ab  e'CRyM  (0, 

Expanding  and  invoking^  the^  orthogonality  of  spherical  harmonics,  and 
performing  the  radial  integration  we  get  [12] 


N'<  1/1  N  +  i.ft' N»L' 

I  I  «!ci:'.‘-*'{/./)(frt) 
(-0  1-0 


N'JL’-L  +  I  +  / 


-  f  )r^"  [« w  +7)p) 

-r.  ^  I  /  (~i)' 

*L(;r^\[«(r-K)P  [a(r-f)r"vj’ 


wliere 


/(2L  +  1)(2L*  +  1)(L  +  M)l(i:'-A/)I\ 

I  (2/^')l(2/V)l(Z.’  +  Af)l(Z.-A/)f  / 


Table  III  shows  the  programming  of  the  overlap  formula.  We  have  set  N'  =  nnp, 
L'  =  hhp,  K  =  kk.  For  clarity,  the  overlap  formula  is  assembled  from  parts  to  get, 

S  =  overlap.  To  compare  with  the  literature  [11],  we  set  =  10,  and  ^  =  2,  and  the 

distance  between  orbits  as  1.4.  Then  we  let  a  =  C’(l-4)  =  14,  b  =  C(l-4)  =  2.8  = 
2(14/10).  Mathematica  can  give  arbitrary  precision  for  arithmetic  operations  if  it 
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is  given  numbers  as  the  ratio  of  integers.  This  applies  also  to  roots  and 
exponentiations.  The  operation  N[overlap,  20]  means  for  Mathematica  to  try  and 
get  20  significant  figures  using  the  overlap  formula.  This  it  succeeds  in  doing  as 
confirmed  by  Bhattacharya  and  Dhabal  [11].  With  the  command  N[overlap,  30]  we 
request  30  digits,  and  get  29. 

A  more  revealing  set  of  parameters  for  this  problem  is  a  =  102/100  and  b  = 
101/100.  Now  cancellation  errors  become  evident.  Thus  N[overlap,  20]  yields 
0.4338568005,  which  is  only  10  digits.  Mathematica  trys  not  to  return  worthless 
digits.  The  loss,  we  suspect,  must  be  due  to  cancellation  errors.  This  becomes 
obvious  if  the  positive  parts  and  the  negative  parts  of  the  overlap  formula  are 
programmed  separately  for  20  digits.  Then 

overlap  =  921614991.7515037098  -  921614991.3176469093  =  .4338568005. 

We  may  obtain  20  digits  using  the  command  N[overlap,  30],  which  gives 

.43385680048834139559. 

Hence,  by  the  simple  expedient  of  changing  the  number  of  digits  requested  we 
may  obtain  arbitrary  precision.  Mathematica  can  easily  deal  with  large  numbers 

of  digits,  so  that  any  physical  system  can  be  computed.  (The  exceptional  case  of  C 

=  ^'  (a  =  b)  can  be  programmed  separately  from  a  simpler  overlap  formula  [13]). 
Conclusion 

Recent  investigators  [11,  14,  15]  working  with  overlap  integrals  have  found 
that  more  than  one  formulation  is  needed  in  order  to  cover  the  range  of  useful 
parameter  values  and  avoid  intolerable  cancellation  errors.  This  is  necesary 
when  using  programmes  that  utilize  a  finite  computer  word  length.  In  this 
paper,  it  has  been  shown  that  by  use  of  a  variable  word  length,  it  is  possible  to 
absorb  cancellation  errors  and  still  work  to  a  specified  accuracy. 

We  have  employed  an  Apple  Macintosh  IIx  computer  using  the  1.2  version  of 
Mathematica.  Benchmark  values  for  overlap  integrals  have  been  achieved  (in  a 
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few  seconds).  Other  multi  center  molecular  integrals  are  under  investigation. 
The  general  usefulness  of  our  methods  for  integral  packages  depends  on  newer 
and  faster  versions  of  computer  algebra  schemes  and  the  employment  of 
mainframes  or  workstations. 
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nn=2;hh=l;mm=0;h=0; 


cpol3momial=Sum[a''  (nn+hh+2h-2pp-2vp-2v-k-kp)  * 

(2pp+2v+2vp+kp)  *  (-1)  (v+qp+p+pp+hh)  *  (-1)  '^hh* 
(2hh-2p)  1  *  (2h-2pp) !  /4'^  (hh+h-p-pp)  /  ( (hh-p) !  * 
p!  *pp  !  *q!  *qp  I  *v!  *vp!  *  (hh+mm-2p-q-v) ! )  * 
(nn-hh+2p+2q+2qp) !  /  ( (h-pp) !  *  (h-mm-2pp-qp-vp) !  * 
kp!  *  (nn-hh+2p+2q+2qp-k-kp) ! ), 

{p,  0,  Floor  [  (hh+mm)  /  2] )  ,  {q,  0,  hh+mm-2p}  , 

{v,  0,  hh+mm-2p-q)  ,  (pp,  0,  Floor  [  (h-mm)  /  2] }  , 

{qp,  0,  h-mm-2pp)  ,  (vp,  0,  h-mm-2pp-qp)  , 

{k,  0,  nn-hh+2p+2q+2qp)  , 

(kp,  0,  nn-hh+2p+2q+2qp-k)  ] 


3  +  3  a  +  2a2  +  +  3r  +  3ar  +  2a2r  +  r^  +  ar2 


cmatrix=CoeftlcientList  [cpol3momial,  (a.r)  ]; 
csscmatrix; 


MatrixForm  [c] 
3  3  1 

3  3  1 

2  2  0 

10  0 


Table  I.  Generation  of  the  zeroth  harmonic  of  the  C  matrix  for  the  2p  orbital 
starting  with  the  C  matrix  polynomial. 


— It — 


alpha=Sujni  [c  [  [i+1,  j+1  ]  ]  *  ( (-1)  j*Exp  [r]  -  Exp  [-r] )  *  a'^i* 
(j*h-l) ,  {i,  0,  nn+hh+h)  ,  (j,  0,  nn+h)  ]  ; 

alpha=(2h+l)  *  (h-mm)  !  /  (2*  (h+mm) !  )  *alpha  ; 

alphae=Series  [alpha,  {r,0,5)  ]  ; 

e=CoefficientList  [alphae,  (r.a)  ]  ; 

MatrixForm  [e] 

0  0  0  1 

0  0  0  0 

0  0  -2/3  1/6 

0  0  0  0 

1/15  1/15  -1/15  1/120 


Table  II.  Generation  of  the  zeroth  harmonic  of  the  E  matrix  for  the  2p 
orbital.  The  alpha  function  polynomial  is  expanded  in  a 
Taylor  series  in  r,  and  the  coefficients  are  collected  in  matrix 
form. 


nnp=l;  hhp=0; 
a=14;  b=2*  (14/10) ; 


kk=2'^  (nnp+nn)  *  (-1)  '^mm*Sqrt  [  (2hh+l)  *  (2hhp+l)  *  (hh+irun) !  * 
(hhp-mm) !  /  ((2nnp)  !  *  (2nn) !  *  (hhp+mm)  I  *  (hh-mm) !  )  ] 

sl=Sum  [  (nnp-hhp+j) !  *c  [  [i+1,  j+1]]  (nnp-2hhp-hh+i+j)  * 

((-1)  /  (a-b)  (nnp-hhp+j+l)  -1  /  (a+b)  (nnp-hhp+j +1)) , 

(i,  0,  nn+hh+hhp)  ,  (j,  0,  nn+hhp)  ]  ; 

s2=Sum  [(nnp-hhp+j)  !  *c  [[i+1,  j+1]]  *b'^  (nnp-2hhp-hh+i+j)  / 
(nnp-hlip+j-k) !  *  ((-1)  i  /  (a+b)  (k+1)  -  (-1) 

{i,  0,  nn+hh+hhp)  ,  (j,  0,  nn+hhp)  ,  (k,  0, 

sl2=Exp  [-b]  *  sl+Exp  [-a]  *s2 

overlap=kk*  (a/b)  (nnp+1/2)  *  sl2  ; 

N  [overlap,  20] 

0.11741378968662828485 

N  [overlap,  30] 

0.11741378968662828485490731401 


Table  III.  The  implementation  of  the  overlap  formula  to  find  S(ls,  2p).  The  Is 
orbital  has  a  screening  constant  of  10,  and  the  2p  orbital  has  a 
screening  constant  of  2.  The  orbitals  have  a  1.4  unit  separation. 
Mathematica  is  requested  to  produce  a  20  digit  result  and  a  30  digit 
result. 
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Abstract 


The  Lowdin  alpha-function  method,  in  which  displaced  orbitals  are 
expanded  in  an  infinite  series  of  spherical  harmonics,  is  implemented  for 
Slater-type  orbitals  using  a  commerical  computer  algebra  program, 
Mathematica.  The  program,  which  is  included,  generates  a  C  matrix  with 
integer  elements  that  characterizes  our  approach  to  multicenter  molecular 
integrals.  The  general  two-center,  two-electron  Coulomb  repulsion  integral  is 
produced  analytically  with  a  finite  number  of  terms.  Each  Coulomb  formula  may 
be  evaluated  to  arbitrary  precision,  since  Mathematica  works  with  integer 
arithmetic.  Hence,  cancellation  errors  can  be  overcome. 
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Introductioii 


The  difficulties  of  working  with  Slater-type  orbitals  (STOs)  when  doing 
variational  treatments  of  molecules  are  proverbial  [1].  Nevertheless,  progress 
continues  to  be  made  [2] .  The  hope  is  that  the  integrals  resulting  from  use  of 
STOs  will  be  efllciently  done  so  that  their  good  physical  characteristics  can  be 
utilized. 

In  this  paper,  two-center,  two-electron  repulsion  integrals  of  the  Coulomb 
type  will  be  evaluated  to  specified  accuracy  by  use  of  a  commerical  computer 
algebra  program  called  Mathematica  [3].  Our  method  is  that  of  the  Lbwdin 
alpha-function  [4],  in  which  displaced  orbitals  are  expanded  in  an  infinite  series 
of  spherical  harmonics  with  functional  coefficients  (alpha-functions).  By  use  of 
an  in-house  version  of  computer  algebra  we  have  been  able  to  characterize  each 
STO  by  a  C  matrix  [5].  The  fact  that  these  matrix  elements  are  integers  is  the  key 
to  our  being  able  to  take  advantage  of  the  high  level  language  and  power  of 
Mathematica. 

We  have  already  shown  [G]  how  Mathematica  can  generate  a  C  matrix  for 
each  displaced  STO  and  how  a  formula  for  overlap  integals  can  be  evaluated  to 
arbitrary  accuracy.  The  more  difficult  Coulomb  problem  involves  the  integration 
over  two  electrons  and  thereby  casts  light  on  how  all  multicenter  molecular 
integrals  may  be  done. 

The  Alpha-function  and  C  matrix 

The  Lbwdin  a-function  method  expands  displaced  orbitals  in  an  infinite 
series  of  spherical  harmonics  about  a  single  center  (the  origin)  before  the  required 
integration  is  attempted  [4].  The  functional  coefficients  of  the  spherical 
harmonics  are  designated  as  a-functions.  We  take  a  displaced  STO  to  be  centered 
at  (0,0, a)  in  its  local  coordinate  system  (R,  0,  (]>),  and  write  it  in  terms  of  the 
original  coorinate  system  (r,  0,  <]))  [5,7].  Thus, 
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;i:  =  AR"^-^e-^*Y^(e,<^). 


a:  = 


'N-l 


■(2L4- IXL  +  M)!' 

2 

4zr(^-f-M)! 

47r(L-M)! 

2!Lj 

t=M 

(-1) 


M 


X  (Ca,Cr)Y“(6>,0). 


where 


(  2«  +  M)!  i^o  j-^0  ' 


X  H..(fa)‘“^“^“^(fr)j  '  ^ 


and 


[(-ly  e^'-e'^']  ,  r<a 
[(-1)‘  e^“-e"^“]  ,  r>a 


A  =  (2^)N+1/2  [(2N)!]'^2  15  (,iie  normalization  factor,  N,  L,  and  M  are  the  quantum 
numbers  of  the  orbital,  and  ^  is  the  screening  constant  or  orbital  exponent. 

Most  importantly  for  our  developments  the  elements  of  the  C  matrix  are 
integers.  Originally,  they  were  obtained  by  programming  the  following 
expression,  using  FORTRAN  and  a  simple  in-house  version  of  computer  algebra 
[5,7]: 
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I(L+M)/2]  L-I-M-2P  L+M-2p-q 

S  S  cr'(ij)a‘r'  =  S  S  S 

1=0  j=0  p=0  q=0  v=0 

[(/-M)/2]  /-M-2p'  ^-M-2p'-q'  t  t-k 

X  E  Z  ZEE 

p’=0  q'=0  v‘=0  k=0  k‘=0 

^ _ a (2L  -  2p)!  (2£  -  2p)! _ 

4L./+P-P  !(L  _  p)!  p!  p'i  qi  qM  v!  v’!  (L  +  M  -  2p  -  q  -  v)! 

^  _ (N  -  L  -f  2p  +  2q  +  2q' )! _ _ 

U  -  p')!  (^  -  M  -  2p’-  q*- v’)!  (N  -  L  +  2p  +  2q  +  2q'-k  -  k‘)! 

where 

x  =  N  +  L  +  2^-  2p'-2v*-2v  -  k  -  k’, 
y  =  2p'+2  V  +  2v'+k'. 

and 

t  =  N-L  +  2p  +  2q  +  2q' 

Coulomb  Integral 

The  definition  [8]  of  the  Coulomb  two-center,  two-electron  repulsion  integral 
is 

J= II L(l)  z,(2)  x^Cl)  dv,  dvj 

We  superimpose  Xa(l)  and  X^(l)  at  the  origin,  and  Xp(2)  and  Xj^(2)  at  a 
distance  of  a  along  the  z  axis  (0,  0,  a).  The  superimposed  orbitals  may  be  merged 
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and  expressed  as  a  linear  sum  of  orbitals  upon  expanding  the  spherical 
harmonics. 


The  Coulomb  integral  may  be  written  as 

J=|  ;rc(2)  li(2)  V(?2)  dv, 


with  the  potential  given  as 


The  Laplace  expression  of  rjj  is 


-I  oo  ( 

f-E  I 

^12  ^=0  /!=-/ 


ire  r‘ 

27+1 


Y>‘ie„(p,)Yr(&2>V2X 


where  r<  is  the  smaller  of  rj,  r2  and  r>  is  the  larger. 

Using  orthogonality  of  spherical  harmonics,  the  potential  may  be 
determined.  In  our  developments,  in  order  to  avoid  having  to  distinguish  between 
real  STOs  and  complex  STOs,  and  to  facilitate  comparisons  with  the  literature, 
only  the  case  of  orbitals  with  magnetic  quantum  M=0  will  be  considered.  This  is  a 
mild  restriction  easily  removed  once  the  chief  problem  of  dealing  with  radial 
integrals  is  worked  out.  Thus  the  potential  becomes  [8] 

V(r„0,)=  k.  X(L.,L„f)  V,  (r,)  P,(0,), 
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where 


«•  ' 

^  £■  dr,  r"  e-”-'  +  r'  J“dr,  r.-e'"'  . 

w  2  ^ 


with  =  j  (6)  Pl^ (0)  P,{0)siii  0d0, representing  the  Gaunt  coefficients,  with  <u  =  +  C» 

n  =  N.  +  +  A  m  =  N,  +  N,  -  / - 1,  k,  =  A,A,  [(2L.  + 1)  (2L,  + 1)]''*  / 2.  Note  that  Y^{e,<p)  = 


P,(0). 


We  perform  the  indicated  integration  to  get 


_ .  .  .  1  n! 

V,(r3)  =  ;;^|^-e 


^  ni 
(n-k)! 


-  iiill 


+  r'e-“*  T  l£ 

^  (m-k)!  a>‘ 

The  electron  charge  density  for  electron  2  becomes 


;t.(2)  X\(1)  =  ^  R?-‘  (L„L„L)  P.O,) 


with  k,  =  A,A,[(2L.  +  1K2L,  + 1)]''"  /  2,  N  =  N.  +  N,  - 1,  f  =  f.  +  f. 
Noting  the  a-Ainction  expansion 

RN.ig-tn.  1  £  (0^) 


and  taking  into  account  orthogonality  =  /)  we  finally  get  [9J 

J  =  ^I  I  ^  (L„L„L)  (L..M)  p dr  V, 
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The  sums  over  L  and  I  are  limited  because  of  the  triangular  rule  for  the  Gaunt  coefficients. 
For  convenience,  we  divide  the  integral  into  three  parts  corresponding  to  the 
three  terms  of  the  potential  [10]. 

Let  ij  +  ij+ij  =|r*  dr  a^.  Thus, 


N+L+/  N+/  I 

i.  =  I  S  C,(i,j) 


i=0  j=0 


CO 


x|(-l)^  J‘’dr  +  (-1)*  j  dr  H  -e  dr 


^  ^  ^0  (n-k)! 


-k  g-(a>-C)r 


-(-1)'  e^"  f~dr 

J  a 


-2<+n-k  g-(a»+f)r 


+  6^“  Tdr  } 

J  O 


is  =  S  X  Z  c,(i, j) 

i  j  (m-k)! 


x{(-l)J  e-^"JJdr  +  (-1)‘  JJdr  r 


J+l+m-k  Q-(o»-C)r 


-e-f"  Tdr  } 

«^o 
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Each  term  of  the  summation  over  I  and  L  is  given  by 


123 


=  (ii  +  i.  (L«.Lb,^)  (L.,Ld.L) 


X 


^L+2/+2^L+/+l  » 


the  (2^+ 1)  factor  having  been  cancelled. 

Basic  integrals 

To  carry  out  the  analytical  evaluation  of  the  Coulomb  integral,  the  following 
basic  formulas  are  used  [llj  (we  have  replaced  the  zero  on  the  integrals  by  the 
inflnitessimal  £). 


gl(n,fl)  =  f  dr  r"  = - - - n  -1; 

gl(n,a)  =  fna-^ne,  n  =  -l 


fl(n,a,b)  =  J  dr  r"  e 


.-ba  n 


n! 


n-t  .  n! 


fl(n,a,b)  =  -(-br"'‘e-'’'’  V ‘  ■  "  ----  - - — -  +  ^ Ei(-ba) 

S  (-n-D!  (-ba)-"-*  (-n-1)! 


+  (-b)  ”  ‘  Y  ■  ^  J  ^  [to|b|  +  r  +  tnel  n  <  -1; 

(n-l)!(-n-t)!  (-n-D!  ^  ''  '  ' 


fl(n,a,b)  =  Ei(-bfl)  -  [f n|b|  +  7  +  Me],  n  =  -L 
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To  obtain  the  e  limit,  e'^e  has  been  expanded  and  all  powers  of  e,  except  e®,  have 
been  dropped  [12].  From  physical  considerations,  all  inverse  powers  of  £  and  f  n  £ 
must  cancel,  if  all  parts  of  the  program  are  considered.  The  other  needed 
integrals  are: 


f2(n,a,b)  =  f  dr  r"  e 

Ja 


e-‘“  ^  n! 


f2(n,a,b)  =  -^  S  ^  (b<.r-  .  nao; 


f2(n.<,,b)  =  (-b)--  e-  Y  tlLzlzM  _J__  - 

(-n-D!  (-ba)-"-‘  (-n-1)! 


Ei(-ba),  n<-l; 


f2(n,a,b)  = -Ei(-ba)  ,  n  = -h 


f3(n,b)  =  dr  r"  e 


=  T#r.  naO; 


f3(n,b)= 


^  (-n  -  1)!  (-n  -  t)! 


(-n-1)! 


n|b|+  y  +  /ne] ,  n  <  1; 


f3(n,b)  =  —  [^njb(+7+ /nej,  n  =  — 1. 


For  the  case  of  Coulomb  integrals  [10,13]  we  may  drop  all  Ei  and  logarithm  terms. 
Also,  Euler's  constant  y  does  not  appear.  Hence,  our  Mathematica  progam  does 
not  include  these  terms. 

Programming  of  the  Coulomb  integral 

Mathematica  is  a  high  level  language  that  has  a  remarkably  close 
correspondence  to  standard  mathematical  notation.  Thus 
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imax  jmax 

Yj  ^(expression)  =  Sum  [expression, {i,iinin, imax,  jjmin,  jmax}] 

i=iinin  j= jmin 

Integrate[f(x),(x,xinin,xmax 

For  added  flexibility  for  iteration  we  may  use  Do-Ioops: 

Do  [Do[equations,{i,imin, imax}, {j,  jmin,  jmax}||. 

It  is  convenient  to  use  a  Which  statement  that  deflnes  a  function  for  various 
parameter  values: 

f[n_,x_|:  =  W}iicli[(est  l,fuiiclion  l,(es(  2, function  2,  True,  function  3). 

In  this  case,  if  test  1  and  test  2  are  unsatifled,  function  3  is  selected. 

After  running  the  program  of  Table  I,  for  the  case  of  all  Is  orbitals  with  equal 

screening  constants  (Cg  =  Cj,  =  =1,  w  =  C  =  2,  typing  "coulomb"  produces 

Roothaan's  [13]  formula  (2p  =  Ca  =  2a).  Typing  a  =  1/100  and  N[coulomb,  20] 
produces  0.62499166683332546.  Typing  N[coulomb,30]  produces 
0.62499166633325460009943731.  Although  we  have  requested  20  digits  and  30  digits 
we  get  17  and  27  digits  computed.  This  is  because  of  cancellation  errors  caused  by 
the  differencing  of  nearly  equal  numbers.  Mathematica  trys  not  to  produce 
worthless  numbers. 

To  conform  to  Mathematica  protocol  and  for  clarity  in  typing  the  following 
notation  is  used: 

Na  =  nna,  =  lilia,  ^  =  alpha;  =  nnb,  Lb  =  hhb,  ^b  =  beta;  M  =  mm;  A^  =  aa, 
Ab  =  ab,  Ac  =  ac,  Aj  =  ad;  L  =  hh,  f  =  h. 
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The  Gaunt  coellicients  are  produced  by  explicitly  integrating  over  three  Legendre 
polynomials. 

Results  for  unequal  screening  constants 

Table  II  shows  the  results  for  several  examples  of  Coulomb  integrals.  These 

examples  deal  with  the  more  difiicult  case  (for  other  methods)  of  ^  ^c+  ^d 

(w  zeta).  In  our  programming  we  are  only  required  to  replace  gl(j-2h+n-k,  a)  by 
fl(j-2h+n-k,  a,  w-zeta]  and  gl|j+l+m-k,  a]  by  fllj+l+m-k,  a,  w-zeta].  The  first 
example  produces  the  12  digits  given  by  the  author  [14]  using  a  Taylor  series 
expansion  of  a  computer  algebra  formula  using  p  and  t  (p  =  (w  +  C)  a/2,  t  =  (w  -  C)/ 
(w  +  C);  p  =  0.02,  t  =  0.01).  Hence,  we  see  that  a  closed  formula  can  effectively 
control  the  word  length  of  the  computer.  For  the  case  with  a  =2,  there  is 
agreement  with  the  seven  decimal  digits  obtained  by  using  expansions  of  the 
alpha-function  [9,15]  in  a  FORTRAN  implementation  for  the  Coulomb  integral. 
Finally,  we  take  two  extreme  examples  using  4f,  3d,  and  Is  orbitals  with 
separations  of  a  =  0.01  and  a  =  100.  With  a  =  0.01  we  need  to  request  40  digits  in 
order  to  get  21.  The  example  with  a  =  100  shows  that  small  values  present  no 
problems. 

Conclusion 

The  examples  shown  in  Table  II  with  wide  ranging  parameters  indicate  that 
Coulomb  integrals  can  be  computed  to  arbitrary  accuracy  for  all  physical  systems. 
It  is  important  in  developing  fast  alternative  methods  to  have  available  completely 
trustworthy  benchmarks  for  all  parameter  ranges.  This  is  assured  because 
Mathemaiica  uses  integ  arithmetic. 

Other  multicenter  molecular  integrals  are  under  study  using  symbolic 
programming. 
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(*  basic  integrals  *) 

gl  [n_,  a_]  :=Which[n>=0,  a''  (n+1)  /  (n+1) ,  n<-l,  a"'  (n+1)  /  (n+1) ,  True,  0] ; 

fl  [n_,  a_,  b_]  :=Which[n>=0,  (-Exp[-b*a]  /h''  (n+1)  *Sum[n!  /  (n-t)  !  *  (a*b) 
(n-t), {t,0,n}l+n!/b^(n+l)), 
n<-l, 

((-b)''(-n-l)*Suin[(-Exp[-b*al*(-n-t-l)  !/(-n-l)  ! 
/(-b*a)-(-n-t)  +  (-n-t-1) !/{-n-l) !/(-n-t) !), 
{t,l,-n-l}]). 

True, 0] ; 

f2  [n_,  a_,  b_]  :=Which[n>=0,  Exp[-a*b]  /b^  (n+1)  *Sum[n!  /  (n-t)  !  *  (b*a)  (n-t) 

(t,0,n}], 
n<-l, 

-  (-b)  '' (-n-1)  *Sum[-Exp[-b*a]  *  (-n-t-1)  !/ (-n-1)  !/ 

(-b*a)  (-n-t)  ,  (t,  1,  -n-l)  1 , 

True, 0] ; 

f  3  [n_,  b_]  :  =Which  [n>=0 ,  n !  /b'^  (n+1) , 
n<-l, 

Sum[(-n-t-l) !/(-n-l) !/(-n-t) ! * (-b) - (-n-l) , {t, 1, -n-l}] , 
True , 0 ] ; 

(*data*) 

nna-1 ; hha=0 ;  alpha=l 
nnb=l;hhb=0;  beta=l 
nnc=l ;  hhc=0 ;  gamma=l 
nnd=l;  hhd=0;  delta=l 
inni=0; 

(♦constants*) 

aa=Sqzt  [  (2*alpha) ''  (2nna+l)  /  (2nna)  !  ]  ; 
ab=Sqrt [ (2*beta) ^ (2nnb+l) / (2nnb) ! ] ; 
ac=Sqrt  [  (2*gainina)  (2nnc+l)  /  (2nnc)  !  ] ; 
ad=Sqrt  [  (2*delta) ''  (2nnd+l)  /  (2nnd)  !  ]  ; 
w=alpha+bet a ; 
zeta=ganuna+delta; 

kl=aa*ab*Sqrt [ (2hha+l) * (2hhb+l) 1/2; 
le2=ac*ad*Sqrt  [  (2hhc+l)  *  (2hhd+l)  ]  /2; 
nn=nnc+nnd-l ; 

coulomb=0; 


Table  I.  Program  using  Mathematica  for  the  Coulomb  integral. 


(*do- loops  for  h  and  hh  *) 

Do [  Do [ 

n=nna+nnb+h; 
m=nna+nnb-h-l ; 

cpolynomial=Sviinta'^  (nn+hh+2h-2pp-2vp-2v-k-kp)  * 

r''  (2pp+2v+2vp+kp)  *  (-1)  (v+qp+p+pp+hh)  *  (-1)  ''hh* 
(2hh-2p) ! * (2h-2pp) ! /4^ (hh+h-p-pp) / ( (hh-p) ! * 
p !  *pp !  *q !  *qp !  *v !  *vp !  *  (hh+inm-2p-q-v)  ! )  * 
(nn-hh+2p+2q+2qp)  !  /  ( (h-pp)  !  *  (h-inm-2pp-qp-vp)  !  * 
kp! * (nn-hh+2p+2q+2qp-k-kp) ! )  , 

{p,  0,  Floor  [  (hh+mm)  /2]  }  ,  {q,  0,  hh+iTun-2p} , 

{v, 0, hh+mm-2p-q} , {pp, 0, Floor [ (h-mm) 12]] , 

{qp,  0,h-nun-2pp},  {vp,  0,  h-inm-2pp-qp)  , 

{k,  0,  nn-hh+2p+2q+2qp) , 

{kp, 0, nn-hh+2p+2q+2qp-k) ] ; 
cinatrix=Coef  ficientList  [cpolynomial,  (a,  r}  ]  ; 
c=cmatrix; 

il=:Suin[n!  /w''  (n+1)  *a''i*zeta''  (i+j)  *c  [  [i+1,  j+1]  ]  * 

( (-l)^j*Exp[-zeta*a] *£1[ j-2h,a, -zeta] 

+  (-1)  ''i*Exp  [zeta*a]  *f2  [  j-2h,  a,  zeta] 

-Expl-zeta*a] *f3 [ j-2h, zeta] ) , 

{i, 0, nn+hh+h] , { j, 0, nn+h) ]  ; 
il=Siinplify  [il] ; 

i2-Sum[n!  /  (n-k)  !  /w''  (k+1)  *zeta^  (i+ j)  *a''i*c  [  [i+1,  j+1]  ]  * 

(-  (-1)  j*Exp[-zeta*a]  *gl  [  j-2h+n-k,  a] 

- (-1) ^i*Exp [zeta*a] *f2 [ j-2h+n-k, a, w+zeta] 

+Exp  t-zeta*a] *f 3 [ j-2h+n-k, w+zeta] ) , 

{i,  0,  nn+hh+h) , { j,0,nn+h), {k,0,n}] ; 
i2=Siinplify  [i2] ; 

i3=Sum[m!  /  (m-k)  !  /w*'  (k+1)  *zeta''  (i+j)  *a''i*c  [  [i+1,  j+1]  ]  * 

( (~1) '' j*Exp  [-zeta*a]  *gl  [  j+l+ra-k,  a] 

+  (-1)  ''i*Exp  [zeta*a]  *f2  [  j+l+m-k,  a,  w+zeta] 

-Exp [-zeta*a] *f3 [j+l+m-k, w+zeta] ) , 

(i, 0, nn+hh+h) , ( j, 0, nn+h) , (k, 0,m) ] ; 
i3=Simplify [i3] ; 
il23=il+i2+i3; 
il23^Simplify [il23] ; 

il23=il23/a^  (hh+h+1)  /zeta''  (hh+2h+2)  /zeta''  (nn-1)  *kl*k2*  (2hh+l)  /2* 
Integrate [LegendreP [hha,  x] *LegendreP [hhb, x] *LegendreP [h, x] , 
(x,-l,l)]* 

Integrate [LegendreP [hhc, x] *LegendreP [hhd, x] *LegendreP [hh, x] , 
(x,-l,l)]; 

coulonib=:couloinb+1123 ; 

, (hh, Abs [hhc-hhd] , hhc+hhd, 2 ) , 

(h, Abs[hha-hhb] ,hha+hhb, 2)  ]] 


Table  I  (continued) 
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ABSTRACT 

Using  the  Lowdin  alpha-function  method  in  which  displaced  orbitals  are  ex¬ 
panded  in  spherical  harmonics,  two-center,  two-electron  repulsion  integrals  of  the 
Coulomb,  hybrid,  and  exchange  type  are  done  analytically  using  Slater-type  orbitals. 
Computer  algebra  and  integer  arithmetic  are  used  to  obtain  analytic  results  and 
avoid  cancellation  errors  by  the  generation  of  rational  matrix  elements  for  C,  E,  and 
F  matrices  that  are  used  to  express  the  a-function.  The  formulas  for  the  integrals  are 
kept  simple  by  reversing  the  order  of  integration  over  each  part  of  a  split  quadrant. 
Only  two  basic  integrals  are  used  which  are  first  efficiently  evaluated  by  using  look-up 
tables  and  then  used  repeatedly. 


I.  INTRODUCTION 


The  need  for  basis  sets  to  be  constructed  of  Slater-type  orbitals  (STOs)  when 
dealing  with  current  problems  of  molecular  reaction  dynamics  and  variants  of  density 
functional  theory  has  been  voiced  by  several  quantum  chemists.  And  this  need 
is  still  felt  in  the  more  traditional  studies  of  bound  state  molecules.^  The  lack  bf 
analytic  procedures  have  mitigated  against  the  use  of  STOs  in  contrast  to  the  almost 
universally  used  Gaussian-type  orbitals  (GTOs). 

This  paper  deals  with  two-center  two-electron  repulsion  integrals  of  the  Coulomb, 
hybrid,  and  exchange  type.  The  Lbwdin  alpha-function  method^  is  used  in  which 
displaced  STOs  are  expanded  in  spherical  harmonics,  as  was  first  done  by  Collidge.^ 
This  method  has  been  augmented  by  computer  algebra  and  by  the  use  of  matrices 
with  integer  elements,'^  to  avoid  the  ever  present  danger  of  computer  cancellation 
errors. 

In  addition  to  presenting  a  method  for  doing  all  two-center  integrals  with  the 
possibility  of  high  speed  and  accuracy,  our  method  has  techniques  that  can  be  used 
in  the  solution  of  three-  and  four-center  integrals. 

Much  of  the  early  work  on  two-center  integrals  was  done  in  elliptical 
coordinates.^"®  Later  work  has  used  a  variety  of  more  genera!  methods.® 


II.  ALPHA-FUNCTION  REPRESENTATIONS 
Every  displaced  STO  may  be  expanded  in  an  infinite  series  of  spherical  harmon¬ 
ics;  the  functional  coefficients  being  designated  as  a-  functions.  Assume  that  a  local 
coordinate  system  (i?,  0,  yj)  is  displaced  a  distance  a  along  the  z-axis.  In  terms  of  the 
original  coordinate  system  {r,9,ip)  we  have^® 
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and 


fi  — iP  e^’’  —  e  r  <  a 

e“^'‘[(  — 1)*  e^“  —  r>a 


(3) 


The  normalization  constant  A  =  (2^)^'*’^^^  [(2iV)!]“^/2 ;  N,  i,  and  M  are  the 
quantum  numbers  of  the  orbital;  ^  is  the  screening  constant  or  the  orbital  exponent. 

For  small  values  of  the  parameters  it  is  necessary  to  expand  the  exponentials  in 
the  a-function,  and  by  use  of  computer  algebra  a  triple  sum  is  reduced  to  a  double  sum 
with  an  appropriate  E  oi  F  matrix  with  rational  elements. Using  these  expansions 
our  method  is  stable  for  arbitrary  values  of  screening  constants.  Thus, 


a/(r)  =  I 


^  N+L+l  JMAX  .  r  ,  , 

1=0  j=l 


(4) 


,  /MAX  , 

E  E  Ft(iJ)(Car(Cry-‘-^, 

i=l  j=Q 


r  >  a 


In  out  examples  JMAX  =  I  MAX  =  36  was  sufficient. 

A  simplified  working  form  of  the  a-function  is  obtained  by  leaving  only  the  r 
variable  intact,  and  storing  one-dimensional  V/(j)  and  Zi{j)  matrices. 
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=  < 


(  JMAX 

E  »■  <  a 

j=l 

e  E  ^  r>a 

j=0 


(5) 


III.  COULOMB  INTEGRAL 

The  definition^  of  the  Coulomb  two-center,  two-electron  repulsion  integral  is 


J  =  yy^Xa(l)  Xt(l)  ^12^  Xc(2)  Xd(2)  dui  du2  (6) 

We  superimpose  Xa(l)  and  Xh(l)  at  the  origin,  and  Xc(2)  and  Xdi^)  ®  distance 
of  a  along  the  z  axis  (o,  o,  a).  The  superimposed  orbitals  Xc(2)  and  x<i(2)  may 
be  merged  and  expressed  as  a  linear  sum  of  orbitals  upon  expanding  the  spherical 
harmonics. 
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The  Coulomb  integral  may  be  written  as 


-  J  >^c(2)  Xrf(2)  V  (f2)  dv2 


with  the  potential  given  as 


V(^2)  =  J  Xa(l)  X6(l)/  ri2  dvi, 


The  Laplace  expansion  of  is 

1  ^  /I  ^ 

-  =  E  E  (9) 

A=0  p=-A  > 

where  r<  is  the  smaller  of  rj,  r2  and  r>  is  the  larger. 

Using  orthogonality  of  spherical  harmonics,  the  potential  may  be  determined. 
In  our  developments,  in  order  to  avoid  having  to  distinguish  between  real  STOs  and 
complex  STOs,  and  to  facilitate  comparisons  with  the  literature,  only  the  case  of 
orbitals  with  magnetic  quantum  A/  =  0  will  be  considered.  This  is  a  mild  restriction 
easily  removed  once  the  chief  problem  of  dealing  with  radial  integrals  is  worked  out. 
Thus  the  potential  becomes® 

^ir2,h)  =  h  ^(La,  Aj,,  A)  V;^  (r2)  Pa(^2).  (10) 

where 


Va(^2) 


^  r  dn  rj‘  dri  rf 


with  (Lo,Lj,,A)  =  /  Pl^{0)  Pl,,{^)  PxW  s\n  0d6,  representing  the  Gaunt  coefRcients, 
and  w  =  Ca  +  Cbi  w  =  Na  +  Nh  +  X,  m  =  Na  +  Ni,  -  X  -  l,ki  =  AaAf,  [{2La  +  1) 
X  (2L6  +  l)l2/2.  Note  that  Y;^{d,^)  =  yJ^PliO). 

We  perform  the  indicated  integrations  using  the  formula 
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+  V  -  — — 
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The  electron  charge  density  for  electron  2  becomes 

x42)  x3(2)  =  I  R?-'  e-<*=  E 

Lf 

with  k2  =  i4c l(2I,c  +  +  1)]^  /2,  iV  =  iVc  +  fVj  —  1,  C  =  Cc  +  Cd- 

Noting  the  a-function  expansion 


rN-I  ,-<R2 


and  taking  into  account  orthogonality  {li  =  A)  we  finally  get 


(14) 


(15) 


j  _  ^1^2 


_  (2L  +  1) 

C^-1  ^  ^  (2A  +  1) 

Ij  a 


i^C}  L)  {La,Lf,,X)  J  r^dr  V;^  (16) 


The  sums  over  L  and  A  are  limited  because  of  the  triangular  rule  for  the  Gaunt 
coefficients. 

The  Coulomb  integral  is  reduced  to  evaluating  basic  integrals  of  the  type 
ra  roo 

I  e~‘^®i”da:,  n  >  0,  and  /  dx,  n  >  0,  or  n  <  0  (17) 

Jo  Ja 

once  the  appropriate  a-function  representation  and  are  substituted.  These  evalu¬ 
ations  are  done  rapidly  using  look-up  tables  and  computer  memory  (see  Appendix). 
To  be  explicit: 
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-A+j-2A  g-(u;+C)r 


ik=0 
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-  1 

Tl  =  iVa  +  +  A 

m  =  iVa  -1-  A^b  —  A  —  1 
=  Ca  +  Cb 
C  =  Cc  +  Cd 


(18) 


(19) 


IV.  HYBRID  INTEGRAL 

All  two-center,  two-electron  repulsion  integrals  are  given  by  the  same  formula 
but  they  differ  in  the  location  of  the  orbitals.  The  hybrid  integral^  hais  three  orbitals, 
Xo(l)>X6(l)  Xc(2),  placed  at  the  origin  and  one  orbital  Xdi^)  placed  at  (o,o,a). 

With  this  being  the  case  the  potential  for  the  hybrid  integral  is  identical  to  that  of 
the  Coulomb  integral.  The  charge  density  Xc(2)  is  somewhat  simpler. 


Xc(2)  x3(2)  =  ^  E  “I.  (’■2)  PU  (»2)  Pl.  («2)  (20) 

Q  li 

Integrating  the  density  over  the  potential  we  get 


I  = 


Aj  A2 


E  E  (i-.ik.-')  (ic.A.U  /  T^dr  V>  aj 

u  X 


(21) 
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Again,  the  sum  is  limited  by  the  triangular  rule  for  Gaunt  coefficients,  and 
upon  a-function  and  V;^  substitutions  the  hybrid  integral  reduces  to  evaluating  basic 
integrals.  To  be  explicit: 


Q  A  U 

fJMAX  ,  fa 

j=U 


+  y  "’(-1)  -  ydj)  r  ^^c+n-t+j-A  g-(Cc+u.)r^^ 

1“  (n  —  ^  JQ 
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r  \  yoo 
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(22) 


Symbols  are  as  in  the  Coulomb  case. 


V.  EXCHANGE  INTEGRAL 

For  the  exchange  integral^  we  place  Xo(l)  and  Xc(2)  at  the  origin,  and  Xi(l)  and 
Xdi"^)  (ojO)®)-  We  substitute  in  the  following  expressions,  invoking  orthogonality, 
and  obtain  an  infinite  series  for  the  exchange  integral,  K. 


Xa{l)  =  Aa 

^6  h  ^  ' 

Xc(2)  =  Ac  rf‘-‘  =-<=’■■  rl(«2.V2) 
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J  j  dr\  dr2  e  af^  e  af^  (r2) 


(24) 


On  a  previous  occasion we  have  first  determined  the  potential,  but  this  leads 
to  di'  cult  basic  integrals.  To  avoid  this  we  shall  split  the  (rj,  r2)  quadrant  as  in 
Figui  I,  as  was  done  by  Lundquist  and  Lowdin.^^  For  arbitrary  functions  /(r^)  and 
g{r2)  we  get  a  more  symmetry  integration  scheme.  Thus, 


noo  1  X 

dri  dr2  /(ri)  g{r2)  -3^  =  /  dr2  p(t’2)-Xh  /  ’'l 

Ty  Jo  T2  Jo 

fOO  1  rr, 

+  /  dri  /  dr2  9{r2)T2 

*/  0  7*  J  V  0 


Therefore,  we  may  write 


(25) 


A-,  =  A-;  +  K{‘. 


(26) 


W1 


ith 


{Lc,X,la) 

%  „  1  r2  (27) 

X  /  dr2  r2‘^^  e  e  aj^(ri)  rf . 

•/  0  ^2 

In  a  analogous  manner  v,  ritten. 

To  make  our  method  analytic,  the  first  integral  for  the  "potential  is  done  by 
a  simple  in-house  version  of  computer  algebra.  Because  two  representations  of  the 
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a-function  must  be  used,  depending  on  whether  r  is  less  than  or  greater  than  a,  six 
regions  of  the  quadrant  must  be  addressed  separately.  The  method  will  be  sufficiently 
illustrated  by  just  considering  the  first  three  regions,  that  is 


K{  =  +  A'J  +  Kl 

Region  1.  Here  we  have  r  <  a.  Hence, 

r 


,  fr, 

=  ^  E  dr, 

^2  j=h 

Using  the  formula® 


fr  °° 

/  dx  x^  e-‘^®  =  r”+l  7i!  V 

•^0  k=Q 


k  k 
O)*  r* 


then. 


JMAX  oo 
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El(n  +  fc  +  1)! 


'■k  „Na-\-l+j  +  k 
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where  n  =  //a  +  1  +  J  +  A. 

For  our  examples,  we  take  the  largest  power  of  r2  to  be  36.  The 
generates  and  stores  the  coefficients  of  r2.  Thus, 

36 

^  Ci(A,m)rJ‘, 

m 

The  “potential”  for  Region  2  is  simply 

,2,..>_C2(A) 


rr2,_  ^ 

Uj(r2)  -  TtiT 
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(30) 

(31) 

computer 
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(34) 
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For  Region  3  we  have 


^'>(>•2)  =  iiVr  r  ’■r'*''  '■i  •  (35) 

T’o  -^a  ■_„ 


Using  the  formula 


with  u>  =  Ca  +  Cb  and  n  =  Na  +  j  +  X  ~  If,  v/e  get 


+  e— 


n!  g"-* 
(n  —  A;)! 


5'>V2)  =  pEE^-4.(2)^^r' 


+  4tEE'“"'‘4.(2) 


and  by  computer  algebra 


'2  j  k 


2Na  +  X 


{n~k)\ 


'  f'l(^2)=^  E  C3„(A,n.)rJ'+%f?i,  (38) 

^2  m=0  ’’2 

where  C2R{X,m.)  and  C'3^(A)  are  stored  arrays  of  coefficients. 

Finally,  substituting  in  the  proper  a-functions  for  the  second  integration  we  get 


36  _  ra 

=  E  E  (j)  r^c  +  l+;+m  g-(C.+<c)r 

;=1,  m  -^0 

kI=  Y,  2frf(i)C'2(A)  /  dr 

j=0 

Nd+li  2Na  +  X  -OQ 

^A=  E  E  Zil(j)C3R{A,m)  /  Jr  e-(“+«^".-l+i-/--A 

j=0  m=0 

^oo 

+  E  /  dr 

.•_n  •^o 


These  integrals  are  readily  identified  as  composed  of  our  two  basic  integrals. 
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VI.  DISCUSSION 


Table  I  shows  the  results  of  evaluating  two-center  Coulomb,  hybrid  and  exchange 
integrals  using  various  combinations  of  Is,  2s,  and  2p  orbitals  with  arbitrary  screening 
constants.  Our  method  is  completely  stable  for  these  permutations,  as  is  to  be 
expected.  The  Coulomb  and  hybrid  integrals  are  assembled  from  a  Unite  number 
of  terms  as  determined  by  the  triangular  rule  for  Gaunt  coefficients;  the  exchange 
integrals  are  approximated  to  7  decimal  digits  by  using  12  harmonics.  These 
calculations  were  made  on  a  CDC  CYBER  850  computer.  E  and  F  matrices  and  look¬ 
up  tables  for  An{z)  and  En{z)  are  considered  as  data  which  is  put  into  fast  memory 
at  run  time.  The  “set-up”  time  or  “overhead”  is  in  the  generation  of  Y  and  Z  one¬ 
dimensional  matrices,  the  production  of  basic  integrals,  and  Gaunt  coefficients.  This 
overhead  requires  about  3  s  of  Central  Processing  Unit  time  (CPU)  for  the  exchange 
integrals  and  about  2  s  for  the  Coulomb  and  hybrid  integrals.  Some  of  the  same  values 
are  needed  for  all  of  these  integrals.  In  a  SCF  calculation,  the  set-up  time  would  only 
be  needed  once  for  a  large  number  of  integrals.  The  CPU  time  of  the  subroutines  for 
the  Coulomb  and  hybrid  integrals  was  about  0.25  and  0.4  s,  respectively.  The  CPU 
time  to  generate  12  harmonics  for  the  exchange  integrals  was  about  0.7  s. 

CONCLUSION 

The  analytical  method  outlined  here  for  two-center  STO  integrals  can  be  imple¬ 
mented  with  parallel  processing  or  vector  processing.  And  with  improved  program¬ 
ming  it  should  be  possible  to  significantly  reduce  computer  time.  In  addition,  analytic 
and  semi-analytic  methods  for  three-  and  four-center  integrals  for  STOs  can  use  some 
of  the  strategies  developed  here. 
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APPENDIX 


The  speed  of  our  method  depends  on  the  rapid  evaluation  of  the  baisic  integrals. 
These  basic  integrals  occur  many  times,  so  generally  speaking,  they  need  only  be 
evaluated  for  the  highest  harmonic  and  then  stored  in  fast  computer  memory  to  be 
used  repeatedly.  The  initial  evaluations  of  the  basic  integrals  are  done  using  look-up 
tables. 

All  the  basic  integrals  needed  can  be  written  in  a  simple  form  after  an  obvious 
change  in  variable.  Thus,^ 

Eni^)  =  f  e~^®a:”£ii,  n>0  (l) 

Jo 

and 

An(2)  =  J  e  i”  dx,  n  >  0  or  n  <  0  (2) 

These  basic  integrals  are  accurately  evaluated  over  a  z-grid  with  Az  =  0.1  over  a 
range  sufficient  to  cover  all  the  electron  repulsion  integrals  to  be  investigated.  These 
tables  are  stored  on  magnetic  tape  and  put  into  fast  memory  at  run  time. 

Silver  and  Rudenberg®  showed  how  look-up  tables  may  be  generated  and  En(z) 
evaluated  by  interpolation  using  a  Taylor  series  whose  derivatives  are  given  by  a  shift 
upward  in  n. 

We  found  that  this  approach  can  also  be  used  for  An{z)  even  when  n  is  negative. 
The  simplicity  of  this  procedure  is  shown  by  first  writing  the  /th  derivative  of  En{z) 
and  An(z). 


Then 


~  En(z)  =  (-1)'  En+l(z) 
^  An(z)  =  (-1)^  A„+/(z) 


and 


En(z  -f  /l)  =  En(z)  +  h  + 


dEnjz)  h\ 

dz^  2!  dz^ 


En{z) 


/! 


(4) 


En{z  +  h)  =  En{z)  -  En^i(z)h  d-  + 


+  (-!)' (6) 
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Similarly, 


A.n{z  +  h)  -  An(z)  -  An+i(z)  h  +  An-^2{2^)-^  + - V  (-1)^  An+i(z)—  (6) 

Five  terms  of  the  Taylor  series  give  sufficient  accuracy.  In  our  examples,  n  varies 
from  -22  to  6,  but  we  extend  n  to  10  so  as  to  use  a  full  Taylor  series  in  every  case. 
Hence,  derivatives  need  not  be  stored. 
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reated  separately  to  analytically 


Table  I.  Two-center,  two-electron  repulsion  integrals.  Comparision;  a)  Reference  14; 
b)  Reference  8;  c)  Reference  15. 
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Semi'Analytical  Method  for  FoiuMDenter  Molecular  hit^pnals 
over  Sflatez^l^pe  Qiintals 


Herbert  W.  Jones 
Physics  Department 
Florida  A&M  University 
Tallahassee,  Florida  32307 

Abstract 

A  strategy  for  the  evaluation  of  four-center  molecular  integrals  over 
Slater-type  orbitals  is  developed  using  the  Ldwdin  alpha-function  approach  in 
which  displaced  orbitals  are  expanded  in  spherical  harmonics.  The  harmonic 
potentials  are  produced  analytically  and  evaluated  along  a  grid.  The  harmonic 
charge  distributions  are  given  an  analytical  formulation  and  are  evaluated  over 
the  common  grid  and  numerical  integrations  are  performed,  for  each  harmonic. 
Using  an  example  with  Is  orbitals,  only  nine  harmonics  are  needed  for  good 
results. 

Computer  algebra  and  integer  arithmetic  are  used  to  generate  C,  E,  and  F 
matrices  that  are  stored  as  part  of  the  data  base.  T  and  X  one-dimensional 
matrices  are  introduced  as  an  aid  in  computation.  The  employment  of  look-up 
tables,  and  vector  and  parallel  processing  promises  to  make  this  method,  which 
can  be  generalized,  practical. 
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Introdiictitm 


It  is  generally  agreed  [1]  that  it  would  be  desirable  to  use  Slater-type  orbitals 
(STOs),  because  of  tbeir  more  physical  nature,  in  ab  initio  quantum  chemistry 
and  moleuclar  physics  as  well  as  the  almost  universally  employed  Gaussian-type 
orbitals  (GTOs).  Advances  toward  this  long  sought  goal  have  recently  been  made 
[2],  but  a  comprehensive  STO  computer  code  competitive  with  GTO  codes  has  yet 
to  be  written. 

The  author  has  continued  to  follow  the  path  of  the  Lowdin  alpha-function 
method  in  which  displaced  orbitals  are  expanded  about  a  single-center  in  an 
infinite  series  of  spherical  harmonics  with  functional  coefficients  (Lowdin 

a-functions)  [3].  This  basic  method  has  been  augmented  by  computer  algebra  and 
integer  arithmetic  to  avoid  the  pitfalls  of  cancellation  errors  due  to  the  finite  word 

length  of  computers  [4j.  A  modified  closed  formula  [5]  for  the  a-function  led  to  a 
C-matrix  with  integer  elements  [6,7].  To  deal  with  small  parameter  values,  the 

exponentials  in  the  a-function  were  expanded  using  computer  algebra  leading  to 
rational  elements  for  E  and  F  matrices  [8]. 

The  a-function  method  has  been  applied  to  four-center  integrals  to  obtain 
formulas  [9]  and  anal3dical  procedures  [10].  However,  there  were  several  difficult 
basic  integrals  that  had  to  be  dealt  with  and  one  that  had  to  be  done  numerically. 
A  pragmatic  view  seems  to  demand  that  we  abandon  a  formula  approach  or  an 
all  anal3rtic  approach  for  this  case  and  adopt  a  semi-analytical  method  in  which 
the  first  integral  (potential)  is  done  analytically  and  the  second  integral  (energy)  is 
done  numerically,  for  each  harmonic.  (An  all  numerical  method  has  been  done 
[11],  but  its  lack  of  speed  confines  it  to  a  checking  role.) 

That  a  semi-anal}rtical  method  can  lead  to  good 'accuracy  was  demonstrated 
by  McLean  [12]  and  Clementi  [13]  using  elliptical  coordinates  for  an  analytic 
potential,  followed  by  two  numerical  integrations  for  linear  molecules.  This 
procei^re  was  proved  effective  for  the  general  case  by  Musso  and  Magnasco  [14], 
in  which  a  three-dimensional  numerical  integration  is  required.  The  advantage 

of  the  a-function  method  to  be  presented  here  is  that  it  only  requires  a 
one-dimensional  numerical  integration  for  each  harmonic  used.  This,  together 
with  the  better  systematics  of  the  method,  promises  to  raise  it  to  the  level  of 
practicality. 
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Our  method  will  be  illustrated  by  an  example  using  Is-orbitals.  The 
generalization  to  higher  angular  momentum  is  straightforward  [15]. 


Alpha-Fuimtaon  Representations 


A  Is  orbital  in  its  coordinate  system  is  given  by  X  =  e*CRA/C  If  it  is 
displaced  a  distance  b  along  the  z  axis  it  may  be  expanded  in  spherical  harmonics 
or  Legendre  polynomials  in  the  original  coordinate  system  as 


X  =  C3^  I  a.P,(cose),  (1) 

where 


1+1 

a.  =  21+  1  L 
^  2  i=0 

t+1 

Z  C,(ij)  Hy  (Cb) i- 1- 1  (Cr) 1 
j=0  ^  J 

(2) 

and 

e*Cb 

eC*"  -  e'Cr],  r<b 

(3) 

Hy  = 

e-Cr  [(-l)i  eCb  -  e'Cb],  r>b. 

(4) 

The  E  matrix  results  when  eC^^  and  e'C^  are  expanded  and  the  triple  sum  is 
reduced  to  a  double  sum  by  computer  algebra.  The  F  matrix  results  when  eC^  and 

e’C^  are  expanded.  In  our 

example,  a  36  term  expansion  was  found  to  be 

sufficient. 

Thus 

1 

1+1 

/e-Cb  Z 
i=0 

36 

Z  E.(ij)  (Cb)i--t+l  (Cr)j,  r<b 

(5) 

ai  = 

r 

e-Cr  Z 
j=o 

36 

Z  (i j)  (Cb)  ^■t+  1  (Cr) J-  -^-1 ,  r >  b 

i=t 

(6) 
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A  further  simplification  results  by  just  keeping  r  intact, 

.  ^+1 

Y.(j)  =  CJ  2  E.GJ)  e-Cb 

1=0 


i=^ 


leading  to 


Y^O)  rj  . 

jrl 


r<  b 


a  #  = 


.  ^+1 

e-Cr  Z  Z,0)  rJ-'t-l, 
j=0 


r>  b 


When  it  is  necessary  to  evaluate  the  a-function  over  a  grid  for  large  values  of 

Cb  it  is  expedient  to  introduce  two  more  one-dimensional  matrices  immediately 
derived  from  the  general  C  matrix: 


T^(j)  =! 


2(/4M)  ! 


N+L+^ 

e-Cb  I  cj“  (ij)  (gb)i*L-M 

i=0 


r<  b 


X^G)  =(2^+  I  c>J^(ij)[(.l)i-e-Cb]  (;b)i-L-^-l  ,  r>b. 

2(^+M)!  i=<) 


Hence, 


£  T, O') [(•!)>  eCr-  e’Cnri-t-l,  r<  b 

j=0  ^ 


N*-l 

e<rz  X/(j)rj-f-l, 
j=0 


r>  b 
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^  These  various  representation’of  the  a-function  considerably  simplify  our 
programming  and  also  speed  up  its  execution. 

FouivCenter  Integrals  for  Is  Orbitals 

The  formalism  for  the  radial  part  of  any  STO  is  essentially  the  same, 
therefore  we  may  sufficiently  illuminate  our  method  by  use  of  an  example  using 
Is  orbitals,  given  by  Trivedi  and  Steinbom  [16].  In  this  case.  Is  orbitals  with 

screening  constants  of  ^  =  1  are  located  as  follows: 

Xa(l)  at  (0,0,0),  Xb(l)  at  (0,0,b),  Xc(2)  at  (c,0,0),  and  Xd(2)  at  (0,d,0),  with  b  = 

c  s  d  =  1.0.  The  method  to  be  shown  is  valid  for  arbitrary  b,  c,  and  d,  as  well  as 

^1,,  and  Cd-  Because  of  our  use  of  the  E  matrix  expansion,  our  methods  are 

stable  for  nearly  equal  or  equal  values  of  Ca  and  Since  our  method  is  partly 

numerical,  no  problems  can  originate  with  the  relative  values  of  and  Cd' 

Our  task  is  to  evaluate  the  integral 

I  ss  ^{2  dvidv2  (15) 

in  which  each  orbital  is  at  one  of  four  separate  locations.  We  orient  the  molecule 
so  that  X  a  at  the  origin  and  X  (^)  is  along  the  z  axis  at  a  distance  of  b  from 
the  origin.  Using  potential,  we  have 

I  =  /  V  Xc(2)*  X,j(2)  dv2  (16) 


with 


V(r2)  =  J  X  a  (1)*  X  b  (1)  r^-^i  dvi .  (17) 

Substituting  in  the  Laplace  expansion  for  l/ri2>  the  expansion  for  X  and 
invoking  orthogonality  we  obtain  [10] 
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V(  T2  ,  Og)  =  2  V-  (12)  P.  (cos  Sg), 
I  ^  ^ 


(18) 


with 


V^Crj)  =  4  (;aCb)="  e'C.r, 


(19) 


where  is  the  larger  of  r^  and  r2 ,  and  r<  is  the  smaller. 


Explicitly, 


21*1 


30  j* 

— i - 2  ■5a>  (j)  Pdrr^*i*2  e<T 

if*l  j-4  ^  ° 


36 


+  r2  Z  YJ  (j)  J  dr  ri  +  i-^  e'C^r 


^+1 

+  r2  Z  Zb  (j;  J  dr  ri*2^  V’ 

j=0  t  " 


,  r2<b 


(20) 


Vi(r2)  =  4  (Ca  Cb)^ 


2£+l 


36  5 

1 - Z  Yb  (j)  /  drr-f+j  +  2  e’C^r 

9-MJt  ^  ® 


jZ+l 
2 


/+!  J=0 


^+1  .r, 

Z  Zb  (j)  J  dr  rJ'  +  i 

■t  D 


r<c.+g. 
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We  will  evaluate  each  of  the  harmonic  potentials  along  a  hundred  point 
grid  with  a  spacing  of  0.1  units.  This  is  accomplished,  after  an  obvious  change  of 
variables,  by  use  of  the  following  basic  formulas  [17,18J 


,1  “ 

Jq  x"  e'^  dx  =  n!  e'^  Z 

(iH-k+1)! 


.CO 

x"  e*2*  dx  = 


rr  ui  o-jc. 

I  - 

e-z^J  Z 


k=0  (n-k) ! 


(22) 


(23) 


e~z^^  dx  = 
xD 


(n-k-D!  (-z)l'-i  -  (-z)n-i 
(n-1)!  (tt-1)! 


Ei(-z) 


Next,  we  must  numerically  integrate  the  potential  over  the  charge  density. 
Using  spherical  coordinates,  take  the  center  of  to  be  at  (c,  F,  y)  and  X^  to  be  at 
(d,A,8).  Then  [9] 


X-  =  2  Z 

m 

Z 

(4re)i/2 

.  0^,  ^  (r.Y)  (02,9/  (25) 

in=0 

M=-in 

anfl 

and 

CD 

n 

Xh  =  2  Z 

z 

(4n)i/2 

_  aj  Y^(A,5)*  Y^  (02,92)  (26) 

n=K) 

V  =-n 

2IH-1 

Substituting  in  V, 

X^,  and 

^d 

into  the  equation  for  I  and  invoking 

orthogonality,  we  get 

I  =  Z 

(27) 

with 

=  J  rZ  dr 

(28) 

and 
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p  (rp  =  S  Z  2  2  aj  i4n)^  Z  Z  Y^cr.y)* 

^  “  “  a&t-l  2IH-1  (2£+i)H  ^ 

X  ^  (A,5)  <  m,  4  U.O  1  n,v>  (29) 

The  angular  brackets  represent  the  integration  of  the  product  of  three 
spherical  harmonics  [19]. 

In  our  example,  the  values  of  ^  ,  m,  and  n  are  taken  from  0  to  8  and  we  have 

r  s  A  s  5  =  90^  and  y  =  0^.  The  numerical  integrations  are  done  by  Simpson's 
Rule,  obtaining  I  =  0.345538,  which  agrees  with  the  four  digits  supplied  by  Trivedi 

and  Steinborn.  As  an  internal  check  of  our  method,  we  set  F  =  A  =  5  =  y  =  0®, 

and  obtain  I  =  0.50703  which  compares  well  with  the  exact  value  for  the  hybrid 

integral  Ih  =  0.50704  [20].  Table  1  lists  the  harmonics  and  their  sums  for  our  two 
runs.  The  Central  Processing  Unit  time  on  a  Cyber  850  Computer  was  11.4  s  in 
each  case.  The  basic  integrals  were  pre-calculated,  stored,  and  reused.  The 
matrices  C,  E,  and  F  are  considered  as  part  of  the  data  base. 

Conclusion 

A  feasible  strategy  has  been  formulated  to  evaluate  the  general  four*center 
molecular  integral  using  Slater-type  orbitals.  Both  objectives  of  sufficient 
accuracy  and  speed  appear  within  reach.  Working  within  the  framework  of  the 

Lbwdin  a-fimction  method,  careful  elimination  of  computer  cancellation  errors  by 
use  of  computer  algebra  has  proven  decisive.  The  rapid  and  accurate 
implementation  of  programmes  is  assured  by  pre-calculated  exact  matrix 
elements  and  look-up  tables  as  well  as  computer  vector  and  parallel  processing. 
Improvements  in  algorithms  and  numerical  integrations  are  under  development. 
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Harmonic 


Exchange  integral 
Orbital  locations 
r  =  A  =  5  =  90°  Y=  QO 


Hydrid  integral 
Orbital  locations 
r  =  A=  y=6=0O 


L 


0 

1 

2 

3 

4 

5 

6 

7 

8 


.3465644 

.0000000 

-.0010363 

.0000000 

.0000108 

.0000000 

-.0000004 

.0000000 


.4573980 

.0444000 

.0045491 

.0005674 

.0000911 

.0000188 


.0000014 

.0000005 


Sum 


.3455385 


.5070311 


Table  I.  Four  Is  orbitals  are  located  by  polar  coordinates  as  indicated: 

Xa(0  at  (0,0,0),  Xb^^)  (1.0,  0,0,),  Xc(2)  at  (1.0,  F,  yX  and 

X  (j(2)  at  (1.0,  A,  5).  Screening  constants  are  equal  to  1.0.  The 
harmonics  and  their  sums  are  given  for  the  two  cases. 
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Multicenter  Molecular  Integrals  Using  Harmonic 
Expansions  of  Slater-Type  Orbitals  and 
Numerical  Integrations 

HERBERT  W.  JONES  and  BABAK  ETEMADI 

Deparmtnt  cf  Physics,  rlorida  AAM  Unirrrsily.  Tallahassee.  Florida  J2JQ7 


Abstract 

Formula  and  analytic  meihods  have  previously  been  explored  for  the  evaluation  ot  multicenter  molecular 
integrals  over  Slater-type  orbitals  by  employing  the  Ldwdin  n-runction  approach.  These  procedures  are 
greatly  simplified  by  numerical  integrations.  The  programming  of  this  numerical  approach  is  straight- 
forward  and  hence  can  serve  as  a  check  on  future  developments. 


Introduction 

As  is  well  known  [1],  (lie  general  problem  of  the  evaluation  of  multicenter  mo¬ 
lecular  integrals  using  Slater-type  orbitals  (STOs)  in  basis  sets  is  still  not  competitive 
with  the  use  of  Gaussian-type  orbitals  (GTOs)  in  quantum  chemistry.  However, 
investigators  (2,3]  continue  to  strive  to  bring  this  about  because  of  their  belief  in 
the  inherent  superiority  of  STO  basis  sets. 

By  use  of  the  Ldwdin  alpha-function  method  (4],  a  formula  that  involves  an 
assembly  of  a  triply  infinite  sum  of  formulae  has  been  produced  for  the  four-center 
( (he  most  centers  needed )  multicenter  integral  using  I  s  orbitals  and  equal  screening 
constants  (S).  An  analytical  version  of  this  formula  has  been  produced  [3]  (hat  is 
much  easier  to  program  because  all  of  (he  parameters  are  "ground  up,"  except  the 
radial  variable,  to  produce  a  simplified  cc-function.  This  paper  shows  how  the  mul¬ 
ticenter  molecular  integral  problem  can  be  further  simplified  by  use  of  numerical 
integration  (hat  leads  to  very  transparent  programming  and  reduced  computer  time. 
However,  here,  we  only  achieve  about  five-decimal  digit  accuracy.  But,  of  course, 
accuracy  and  time  is  dependent  on  the  numerical  integration  scheme.  Nevertheless, 
for  the  efficient  writing  of  new  programs,  a  first  run  using  simplified  and  reasonably 
fast  programs  is  greatly  desired;  certainly  the  method  to  be  presented  meets  (his 
criteria. 


Expansion  of  Displaced  Orbitals 

To  illustrate  our  procedure,  it  is  sufficient  to  use  all  Is  orbitals,  with  screening 
constants  equal  to  I  (t  =  I  )• 

We  shift  an  orbital  X  =  from  the  origin  to  a  distance  a  along  the  2- 

axis.  Now,  with  respect  to  the  origin  we  have  (6] 
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X  =(»)-•'*  2  a//’/(C0Si>),  (I) 

1-0 

where 

«/  =  ^  Z  2  C,(|.  (2) 

*  1-0  >'0 

and 

e  -K-OV- e-'J.  r<a 

(3) 

e -'J,  f>a 

For  accurate  evaluation  of  tlie  a-function  at  small  radial  distances,  we  must  expand 
exp(r)  and  exp(-r)  for  r  <  a  and  thereby  obtain  a  power  series  representation 
with  coefficients  given  by  an  E  matrix  (7|.  For  the  case  oFr  >  a,  an  expansion  of 
exp(a)  and  exp(-n)  leads  to  an  F  matrix  (Tors-orbitals  the  F  matrix  is  the  transpose 
of  the  E  matrix). 

Thus 

-a  In  JMAX 

-TTT  2)  2  E,{i.j)a*r^. 

^  1-0  i-t 

IMAX  /♦! 

-JTi  2  llW,j)a‘r>, 

^  i-l  i-o 

A  further  simplification  results  by  just  keeping  r  intact.  Defining 

1 

F,(y)=  Z  e-  (5) 

1-0 

and 

IMAX 

Z,{j)=  Z  r,(i.j)a'  (6) 

i-t 

we  may  finally  write 

JMAX 

Z  F,(y)r^ 

j-i 

It  I 

e-'  Z  . 

J-o 

In  the  case  n  =  2,  we  get  convergence  over  the  grid  from  0  to  2  (r  <  a),  and 
from  2  to  10  (r  >  a).  The  grid  values  of  a;  are  accurate  to  12  decimal  digits  by 
taking  JMAX  and  IMAX  to  be  36.  We  choose  a  grid  of  0.1  throughout. 


r  <a 

(7) 

r>  a 


r  <  a 

(4) 

r>  a 
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Examples  of  Various  MuKicenler  Molecular  integrals 

Overlap 

1  lie  simplest  multicciiter  integral  is  the  overlap 

S  =  f  Xa(l)Xi(l)r/v,  (8) 

Locating  Xa  at  the  origin  and  X/,  at  (0,0, n)t  ^ve  have 

=  (9) 

and 

X*  =  (ir)-''^  2  tr/Afeosd)  (10) 

/ 

Using  orthogonality  of  Legendre  polynomials, 

S  =  4  j  drr^e-'a,  (II) 

The  exact  formula  in  this  case  is 

J  =  e'"(M  n  +  flV3)  (12) 

For  fl  =  2,  5  =  0.58645289. 

Using  a  lUO  point  grid  and  Simpson's  rule  for  numerical  integration,  we  get 

Jf*IO 

r/rf^e"rt/  =  0.58645185  (13) 

0 

The  grid  spacing,  0. 1,  and  grid  length,  10,  have  been  chosen  to  achieve  five-decimal 
digit  accuracy,  thereby  including  99.9998%  of  the  charge  with  the  requisite  accuracy. 
The  overlap  is  done  in  order  to  choose  a  proper  grid.  1'his  task  was  accomplished 
with  0.65  s  of  central  processing  unit  (CPU)  time,  on  a  CDC  Cyber  850  computer. 

Tlircc-Ccntci  Nuclenr-Alfraction  Integral  (Lleetroslatlc  Potential) 

These  two  integrals  (three-center  nuclear  attraction  and  electrostatic  potential) 
differ  by  only  a  constant.  Working  with  the  potential,  we  seek  its  value  at  the  point 
(r],  di)  due  to  a  charge  density  given  as  the  product  of  two  orbitals,  X.  located  at 
the  origin  and  X»  at  (0,0,<i): 

no,d2)  =  /rf»iX.(l)X»(l)/r„  (I4) 

We  substitute  the  Laplace  expansion  for  l/ru, 

-  =  4.r2  Z  (2X+ I)-' 

A-Onf-X 


(I5) 
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where  r>  is  the  larger  of  r|  and  and  is  the  smaller.  Recognizing  the  orthogonality 
of  the  spherical  harmonics,  we  get 

d,)  =  2  »'/(r,)/’,(cos  d,)  (16) 

I 

with 


Assuming  r2  =  0.5  and  dj  =  0  with  a  =  2,  we  have,  dropping  unnecessary  subscripts. 


f'liri)  = 


4  I 
2/+  I  f'j“ 


J^O.5 

dr  r^e~'r'ai(r)  + 

0 


JrL 

21+  I 


«/(/•) 


(18) 


Hence, 

i 

r(0.5,  0)  =  2  =  0.45038  (19) 

1-0 


This  value  is  obtained  by  using  six  harmonics  and  numerical  integration.  This  is 
to  be  compared  with  the  exact  value  (8)  of  0.44996. 

The  cru  time  used  was  0,69  s. 

Exchange  Integral 

For  the  two-center  exchange  integral  we  use  the  same  orbitals  and  their  locations 
as  before.  But  now  the  potential  is  needed  at  100  points  along  the  gnd  that  is 
similar  to  the  r,  grid.  Hence, 

■  21^  I '  <20» 

This  requires  compulation  at  10,000  grid  points  formed  by  a  two-dimensional  (r|, 
ri)  lattice.  This  was  done  using  4.62  s  of  cru  time. 

The  exchange  integral  is  given  by 

JJ  Xa(  I  )Xft(  I  )r-%{  1  )XA  I )  dv,  dvj  (21 ) 

where  x,  and  X,  coincide  as  well  as  X»  and  X4. 

The  exchange  integral  is  equivalent  to 

j  dvjXA2)XA2)y{r2,rfj)  (22) 

Using  the  a-function  expansion  of  X»  and  Xj  and  the  Laplace  expansion  of  I /r,i 
together  with  the  orthogonal  properties  of  spherical  harmonics,  one  obtains  (9] 


(23) 
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with 


=  drr^y,e-'a,  (24) 

Using  five  harmonics  and  numerical  integration,  we  get  K  =  0.184190,  which  re¬ 
quired  4.43  s  of  CPU  time.  The  exact  value  obtained  by  using  the  closed  formula 
ofSugiura  (10)  is  AT  =  0.184156. 

Four-Center  Integral 

We  finally  solve  a  four-center  two-eleciron-repulsion  problem  having  three  Is 
orbitals  located  on  a  sphere  of  radius  I,  and  one  at  the  origin  (11). 

T  his  problem  has  been  done  by  formula  ( 5 )  and  also  analytically  (9).  The  orbital 
Xa(  I )  is  located  at  the  origin  and  I )  at  a  distance  of  I  along  the  z-axis.  The 
orbitals  are  used  to  produce  the  potential  as  shown.  Orbitals  X,(2)  and  Xrf(2)  are 
olT-axis  with  center  locations  given  by  spherical  coordinates  (a,y,r)  and  (<i,j,A), 
respectively,  with  a=  i.  Using  the  Legendre  addition  theorem,  the  expansions  take 
the  form 

“  "*  (4ir)'/* 

Xc  =  2  2  2  (25) 

m-0  »i--in  * 

and 

Xd=  2  2  2  Vt— ««n(A,«)*n(i),,»?,)  (26) 

Making  the  proper  substitution,  we  may  write  the  resulting  integral  as  the  sum  of 
the  product  of  radial  functions  hm  multiplied  by  angular  functions  Ai„„.  Thus 

/  =  Z  2  Z  (27) 

i  m  n 

where 

limit  ~  J"  dr  r  ViOniOii  (28) 

and 


A,„„  =  ^  -77,  Z  Z  r%{l\  T)*l';(A.a)<;;i,  ^|/,  0|/i,  r>  (29) 

where  the  angular  brackets  represent  Gaunt  coelficients.  (P,  y)  and  (A,  d)  are  the 
angular  locations  of  the  center  of  the  orbitals.  The  evaluation  of  Ii^  quite  similar 
to  that  of  Ki,  except  that  the  formula  is  not  closed.  Reasonable  answers  are  obtained 
by  taking  /  from  0  to  4  and  m  and  n  from  0  to  6  with  m  +  ii  £  6. 
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To  dclerniine  the  accuracy  of  lliis  integral,  we  set  -y  =  P  =  4  =  A  =  0  and  get  a 
hybrid  integral  whose  value  is  known  (/  =  0.S0704S).  The  numerical  result  here 
is  /  =  0.506284. 


Conclusion 

It  has  been  demonstrated  that  the  multicenter  integral  problem  is  readily  pro¬ 
grammed  and  evaluated  by  numerical  methods.  A  dilfcrent  numerical  integration 
algorithm  might  conceivably  lead  to  acceptable  accuracy  and  increased  speed.  In 
any  event,  the  ease  in  implementation  of  the  herein  described  method,  at  least  as 
far  as  Is  orbitals  are  concerned,  ensures  that  it  can  serve  as  a  first  check  on  any 
new  developments  in  molecular  integrals. 
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Anal3rtical  method  for  two-  and  three-center  molecular  integrals  over  Slater-type 
s-orbitals  using  expansions  in  spherical  harmonics 


Herbert  W.  Jones 
Department  of  Physics 
Florida  A&M  University 
Tallahassee,  FL  32307 

ABSTRACT 

Using  the  Lowdin  alpha- function  method  in  which  displaced  orbitals  are 
expanded  in  spherical  harmonics,  two-center  and  three-center  electron-repulsion 
integrals  of  the  exchange  and  Coulomb  t3T)e  over  Slater-type  s-orbitals  are 
evaluated.  By  means  of  computer  algebra,  analytical  procedures  are 
implemented  and  no  numerical  integration  is  needed.  The  formulas  for  the 
integrals  are  kept  simple  by  reversing  the  order  of  integration  over  each  part  of  a 
split  quadrant.  Orbitals  of  the  Is  and  2s  type  are  used  to  illustrate  the  method, 
which  can  be  generalized. 
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I.  Introductioii 

Almost  all  recent  calculations  in  ab  initio  molecular  physics  and 
quantum  chemistry  are  done  using  basis  sets  made  up  of  Gaussian-t3rpe  orbitals 
(GTOs),  characterized  by  exp(-R2),  because  of  the  ease  of  computation  of  their 
multicenter  molecular  integrals.  It  is  recognized^  that  basis  sets  composed  of 
Slater-type  orbitals  (STOs),  characterized  by  exp(-R),  are  more  representative  of 
physical  orbitals,  however,  their  molecular  integrals  have  the  reputation  of  being 
"intractable".  Yet,  progress  continues  to  be  made  toward  the  solution  of  the  STO 
integral  problem.2 

In  this  paper,  t;he  method  of  the  Lbwdin  alpha-function^  is  used  again^  in 
which  displaced  orbitals  are  expanded  about  a  common  origin  or  single-center  in 
terms  of  spherical  harmonics  with  functional  coefficients  (alpha-functions).  This 
time,  the  analytic  method  is  simplified  by  defining  a  "potential"  in  various  regions 
of  a  split  (ri,r2)  quadrant.  Now  these  "potentials"  do  not  conttnn  the  Ei  function, 
but  only  powers  and  exponentials.  Hence,  the  second  integration  for  energy 
(Coulomb  or  exchange)  can  always  be  done  analytically  and  only  three  types  of 
integration  formulas  are  needed.  As  an  illustration  of  this  method,  two-  and 
three-center  exchange  and  Coulomb-type  integrals  are  done  for  Is  and  2s  orbitals. 

In  the  Lbwdin  a-function  method  all  terms  separate  into  radial  and  angular 
parts;  our  C,  E,  and  F  matrices  for  higher  quantum  numbers  only  differ  in  their 
numerical  elements.  Hence,  s-orbitals  can  serve  as  an  adequate  prototype  for  the 
general  two-  and  three-center  problem. 

A  simple  form  of  computer  algebra  is  used  in  our  approach  so  that  the 
radial  variable  "r"  remains  identifiabie,  thereby  making  the  method  analytic. 
Integer  arithmetic  is  used  to  generate  our  C,  E,  and  F  matrices  so  as  to  avoid 
cancellation  errors  (these  matrices  are  considered  as  input  data). 
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II.  Alpha*Function  Representations 

Every  displaced  STO  may  be  expanded  in  an  infinite  series  of  spherical 

harmonics;  the  functional  coefficients  being  designed  as  a-functions.  For  a 

s-orbital  in  its  own  coordinate  system  (R,0,(p),  X  =  A  e*CR  YO(0,(p),  where  the 

normalization  constant  A  =  (20^"^^  [2N)!]-^2  ^  ^  is  the  screening  constant  or  orbital 

exponent,  and  the  spherical  harmonic  YO(0,(p)  equals  If  the  orbital  is 

placed  at  a  distance  of  a  along  the  z-axis  of  a  (r,0,(p)  coordinate  system  its 
representation  in  this  coordinate  system  is^ 

X  =  A/(CN-1)  I  (2  £+l)-i/2  a  „  (r)  Yo  (0.9) 

£=0  I  ^ 


where 


cl^t)  = 


and 


21±1  2  E  C.(ij)Hij(Ca)i-^*l(Cr)i-'^-l 

2  i=0  j=0 


e-Ca  eCr  -e'Cr] ,  r<a 


e-Cr  [(.l)i  eCa  -e'Ca  ] ,  r>a. 


(1) 


If  the  exponentials  in  the  a-function  are  expanded,  we  get  a  representation  in  E 
and  F  matrices:® 


i 

a^Cr)  = 


„  N4-^  JMAX 

,  e-Ca  Z  Z  E.(ij)  (Cr)J,  r<a 

i=0  j=^ 

IMAX  N+£  ... 

e-Cr  Z  Z  F.  (ij)  r>a. 

i.=  t  j=0 


(2) 
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A  final  simplification  results  by  substituting  in  the  numerical  values  of  the 
parameters  to  obtain  Y^(j)  and  Z^(j),  which  are  coefficients  of  the  polynomials  in  r. 

JMAX 

/  2  Y»(j)rj  ,  T<a 

|j=^ 

NU 

e-Cr  Z  Zj>(j)rj-l-l  .  r>a  (3) 

J=0 

In  this  paper  we  have  used  36  for  EMAX  and  JMAX. 

III.  Three-Center  Exchange  Integral 

The  two-electron  repulsion  integral  is  defined  as 

K  =  J/dvidv2  X^(l)  X^d)  X^(2)  X^(2).  (4) 

The  locations  of  the  orbitals  for  the  three-center  exchange  integral  is  given  in 

Figiu-e  1.  We  let  a  be  the  distance  to  X  ^d)  and  b  be  the  distance  to  Xjj(2);  9  is  the 
angle  between  a  and  b.  In  their  local  coordinate  systems  the  orbitals  are  given  as 


X^(V=  Aa  yj  (ej,  (pj)  %1,(1)  =  Ab  e'PRi  YO  (©i,  <pi) 

Xc(2)=  Ace-Fj  Y0(92,<P2)  Xj(2)  =  Aj  6-5^2  YJ  (©2,  q)2).  (5) 

We  write  the  expansion  of  X^jd); 

Xbd)=AbI  (1)  YO  (0i,9p.  (6) 
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By  use  of  the  Legendre  addition  theorem  we  have”^ 


X^(2)  =  Ad  i  2  (47t)i/2  (2f  +l)-i  ad  (2)  W  (0i,  (p^)  W*  (02,  (?£)  (7) 

£  =0  M= 


Now, 


l/ri2  =  47t  2  2  (2£+l)-l. 

1=0  np  -£ 


ym*  (01.  (Pi)  Ym(e2,(p2), 


where  r<  is  the  lesser  of  rx  and  r2,  and  r>  is  the  greater.  We  will  use 


yo  (0,  (p)  =  ((2  £+l)/47r)i^2  (cos  0). 


Substituting  these  expressions  into  the  formula  for  the  exchange  integral  and 
using  the  orthogonal  properties  of  spherical  harmonics,  we  get 


K  =  2  (COS0) 


where 


K^=  A^//dridr2  r2  r2  e-%  e-Y^z  «?(!)  (2)  f;  ^  (11) 


A^  —  Aq  Ajj  AcAd  (2£  +  1)'2. 
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To  simplify  the  integrations  necessary  for  we  divide  the  (rj,  r2)  quadrant 
into  nine  regions  as  shown  in  Figure  2. 

Now  we  may  write® 


I  II 


as) 


rN^+l  dri  e-cu*!  (1)  rf  , 

(14) 


II 

K. 


=  A,  r 


rN_  +  l 


1 


a 


dri 


a|(l). 


£+1 


C‘ 


a 


‘^(2)rf. 


(15) 


We  take  the  "potential"  in  the  various  regions  to  be 


f  2  dri  P(i)  rf 
0  £ 


dri^d)  rf 


dr,  ^(1)  rf 
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where 


P^d)  =  rN^^l  e-ar^  (1) 

and 


p(2)  =  rNg+l  e-YT,  (2).  (16) 

I  ^  t 


VVe  take  note  of  the  various  expressions  for  the  a-function  in  the  various  regions: 


JMAX 

0^(1)=  Z  Yb(j)rJ  ,  ri<a  (17) 

^  3=1  ^ 
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1 

.  a 

1 

2 

,b 

p^(2) 

2 

>0 

dr2 

p(2) 

= 

la 

dr2 

3 

b 

3 

4 

p^(2) 

4 

K£  = 

A, 

la 

dr2 

pf 

\ 

lb 

dr2 

It 

<; 

dr2 

5 

6 

= 

C 

dri 

p^d) 

6 

8 


K  .  =  A  .  /  dr,  pi>  V 
^  ^  It 


-1’^*  Pf'"  t 


K^=  A^j;  dr,p»  Y 


Three  basic  formulas  for  integrals  are  needed  (the  first  one  given  by 


Silverstone^); 


f  1  ,  ”  (-Z) 

J  X"  e'“  dx  =  Z)  - 

o  k=0 


w  un 

J  x"  e'“  dx  =  e'^  Z 


n  >  0 


k=0  (n-k) ! 


/  e  ™  dx  =  e-Z  Z  (n-k-l){  Ei(-z),  n>o. 

xn  (n-1)!  (n-1)! 


Taking  note  of  the  different  expressions  for  the  a-function  in  the  various 
regions  and  making  the  appropriate  change  of  variable  in  the  basic  integrals,  we 
get  analytic  expressions  for  the  potentials: 


= 


JMAX  <» 

z  z 

3^1  k=0 


k!  (n+k+l) 


Yb  (j)  rN^+j+k+1 


9 


eo> 


JMAX 
e-«ro  I 


iii=^^+£+l 


•On  2 


JMAX 

-2  I  (-a)»^ 

j=^  k=0  k!(n+k+l) 


(-a)>^  Y^(j)  a^a 


with  n  = 


Vp  = 


,£+l  j=0 


Zb(j)  J^2  dri  rn  e-(cc  +  fi)r^ 


n 

=  _  ei(«  +  P)r2  £  2  Zb  (j)  _n! 


j=0  k=0 


(n"*ls) !  ir+i 


n 

2  I  Zb(j) 
j=0  k=0 


n! 

(n-k)! 


a  n-k  Q-(a  +  P)  a 


3  3R  3A 

v^=  V,  ^  V, 


=  e'(«  +  P)  r,  Z 


+  p. 


with  n  =  Nq  +  j. 
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JMAX  a, 

Z  I 
q=£  k=0 


Y^(q)  rN.+q+k+i 


k!  (n+k+1) 


6  JMAX 

V  =  e-TT,  Z 
^  ni=Nc 


JMAX  6 

Z  P 


JMAX 

-Z 

q=l 


i~y)^ 

k!  (n+k+1) 


e-^  Yd(q)  bN  +-t+k+q+2 


with  n  =  Nr  +  +q+l. 


Zd(q)  J  1  drg  rNr'*'*! 

-CD  2  ^ 


n 

e-(l'+S)  Tj  j;  £ 

-  q=0  k=0 

^+1 

1 


(q)  n! 
(n-k)! 


(Y^6) 


J^l 

I  Z 
q=0  k=0 


(q)  n! 
(n-k)! 


bn-k  e‘(Y+5)b 

(Y+6)*^+1 
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9  ® 

V  =  V 
I  I 


SB 

V. 


v^= 


t+2 

e-(Y+5)rj  £ 

1 


SR 

P 


bn 


fin 


9B 


.^+l 


(27) 


with  n  =  Nc+  q 


The  coefiBcients  P 


bn 


2 

P 


3R 


3A 

P« 


6 

P 


tm 


8 

P 


9R 

P 


-fin 


,  and 


3B 


are  generated  by  computer  and  stored. 


Now,  we  put  together  the  final  integrals: 


^  i  a 

K  =  A  I  I  (q)  P  /  dr  e'7**  r  N  +  q+m+1 


=  A^S  YJ  (q)  dj.  e-F  rN,,+  q--C 

q  (L 


3  3R  3A 

0 

3®  3H 

K  =  A  2  ilYMq)  P  /  dr  e'(tt  +  p  +  y)r  rN+q+m-^ 


3A  3A 

=  A^I  Yd  (q)  P  dr  e’F  rVq-^ 
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^  Z  Z  Zd  (q)  J  dr  e-(a  +  p  +  Y+5)r  j-N^ 


»  =  A-Z  Z  YJ>  (j)  f  dr  e-ar  rN^+j  +  m+l 
^  m  t  Ira  o  “ 


K,= 

A 

A.  ^ 

5 

5R 

Kf  = 

5R 

= 

A.  Z 
q 

5A 

K^  = 

A„Z 

q 

6 

= 

7 

= 

A  I 

1 

8 

K  = 

A  Z 

1 

^  j 

9 

9R 

K,  = 

K 

1 

1 

9R 

= 

98 

A  z  : 

j 

^4^  Jb 

5A 

3R 


3A  „ 

^^(q)  P,  1, 


=  A  Z  Z  Zb  (j)  P  dr  e*(a  +  P)r  rN^+j  +  m 


8 

S  Zb(j)  P  f  dr  e-(a+P)r  rN.-l+j- 
I  j  I  I  b  “ 


2^ 


9B 


gR 


;.  =  A  -  Z  Z  Zb  (j)  P  L  dr  e-(a  +  P+y+5)r  j-N 

^  Z.  A  wr»  y  Pm  H  C 


■*-  J  m  b 

3B 


l+q+m-2 


-  £ 


-1  +  j+m-2^ 


(28) 


In  the  summations,  the  maximum  power  of  r  used  was  36. 

As  an  example  of  the  three-center  exchange  integral  we  take  the  case  of  a  = 

b  =  2,  0  =  45®,  and  a  =  P  =  Y=5=  1.2.  Table  I  shows  the  results  of  our 


computations.  The  last  columns  give  the  sum  of  the  preceding  columns.  We 

multiply  each  last  column  value  by  P  (cos  0)  and  add.  Hence,  K  =  .1099365  and  K 
=  .0845324,  and  therefore  K  =  .1944689.  (The  Central  Processing  Unit  time  was  93  s 
on  a  Cyber  760).  This  value  agrees  with  the  six  figures  given  by  Trivedi  and 
Steinbom^O,  and  Graovac,  et  al.^^ 

We  may  readily  specify  the  three-center  exchange  integral  to  the  two-center 

exchange  integral.  As  an  example  we  set  a  =  b  =2,  0  =  0°,  and  keep  all  screening 
constants  the  same.  Because  of  symmetry  considerations  we  may  write  K  =  2K, 
and  the  problem  reduces  to  integration  over  three  regions.  The  value  computed 
was  K  =  .1433972  (17  s  CPU  time).  This  compares  well  with  the  exact  value  of  K  = 
.1433970  obtained  by  use  of  Sugiura’s  formula.^2  This  kind  of  comparison  is  an 
excellent  check  on  our  method. 

Let  us  consider  two  additional  examples.  For  a  two-center  exchange 
integral  take  (1)  to  be  a  2s  orbital  with  Xj,  (1) ,  X^  (2),  and  Xjj(2)  to  be  Is  orbitals, 

with  all  screening  constants  equal  to  1.0.  Let  the  displacement  distance  be  a  =  1.0. 
Then  we  get  K  =  0.35678082,  which  agrees  with  the  seven  digits  given  by  Maslen 

and  Trefry.l3  If  we  interchange  X^  (1)  and  X^^  (1),  we  get  the  same  answer,  which 
demonstrates  the  consistency  of  our  method.  In  the  second  example,  we  study  a 
three-center  exchange  integral  with  all  2s  orbitals,  a  =V%r  b  =  2.0,  0  =  60°,  a  =  0.5, 
P  =  1.0,  Y=  1.2,  and  5  =  1.5.  We  find  that  K  =  0.13960560. 

IV.  Three-Center  Coulomb  Integwd 

Coulomb  integrals  are  simpler  to  formulate  than  the  corresponding 
exchange  integrals.  This  is  because  we  may  merge  the  two  orbitals  of  the  same 
electron  and  place  it  at  the  origin.  Then,  the  potential  due  to  these  orbitals  can  be 
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immediately  written  as  an  analytical  expression.  14.15 

For  the  three-center  Coulomb  integral  we  use  the  same  geometry  as  before 
(Fig.  1),  but  now  we  interchange  the  positions  of  orbitals  (1)  and  X^(2).  In  this 

case,  we  may  consider  the  product  of  the  two  orbitals  at  the  origin  as  a  single 
charge  distribution.  Thus 

X^(l)Xb(l)  =  A3AbrN^+Njj-2e-(a  +  P)ri  YO  (0i,(pi)  /  (47c)i^  (29) 

The  true  potential  is 

V(r2)  =  /dvi  X^(l)  Xj,(l)/ri2  (30) 

The  integration  is  easily  carried  out  getting 


V(r2) 


8  (g  p)3^ 
(ori-B)’^^ 


_  e-(a  +  P)r2 


(31) 


We  write  each  term  as  a  separate  potential 


2 

V(r2)  =  ^-(a  +  P)r2 
^2 

3 

V  (r2)  =  -  (  a  +  p  )  e  +  P)r2 
2 


(32) 


(33) 


(34) 
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The  definition  of  the  Coulomb  integral  is 


J  =  /J  dvidv2  x^l)  r^i  X*(2)X^(,2).  (35) 

In  our  case 

J  =  /  dv2  V  (r2)  \  (2)  X^  (2).  (36) 

Substituting  in  the  appropriate  a-functions  and  taking  orthogonality  into 
account  we  get 


J  =  S  (cos0) 

with 

Jl  =  J  r  2  dv  V  (r)  aj  (r)  oj  (r).  (37) 

This  time, 


Now,  we  simply  write 

12  3 

=  /r2  dr  a®  {V  +  V  +  V  }  (38) 
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Taking  into  account  the  range  of  validity  of  each  a-function,  we  write 


i 


Furthermore,  let  each  region  I,  II,  and  III  define  three 
to  the  three  potentials. 

9  i 


Thus,  finally  J .  =  Z  J 
^  i=l 


(39) 

values  corresponding 


(40) 


Table  II  tabulates  J ^  for  13  harmonics.  The  last  column  corresponds  to  . 

Mtilti plying  the  last  column  by  P^(cos9)  and  adding  we  get  J  =  .343848  (CPU  time 
is  164  s). 

We  may  check  our  procedures  by  setting  o  =  b  =  2  and  0  =  0°.  This  gives  the 
result  J  =  .454950.  Closed  formulas  are  availablel®  for  two-center  Coulomb 
integrals.  For  this  example,  the  formula  gives  J  =  .455049. 
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V.  Two-Center  Coulomb  Integral 


The  two-center  coulomb  integral  forms  an  exceptionally  simple  case  with 
our  method.  With  s-orbitals  only  one  harmonic  is  needed,  namely,  t  =  0  .  This 

comes  about  if  we  merge  orbitals  (2)  and  Xj  (2).  Thus, 


X^{2)X^(2)  =  (47t)-i^2  A,  Ad  e-(Y+5)R2  yo  (0.()))  (41) 


We  consider  this  charge  density  as  essentially  a  Is  orbital  with  a  screening 
constant  of  y+5.  We  write  its  a-function  as  .  Hence, 


X,(2)  Xd(2)  =  YJA,  Ad  I  (2  1+  af  (2)  Yo  (42) 

Substituting  this  expression  into  the  formula  for  the  Coulomb  integral  and 
invoking  orthogonality,  only  the  X  =  0  term  survives.  Thus  we  get 


J  =  YCD(j)j^  drrj  +  2  v+  A„I  Z^D  (j)  dr  rJ  +  i  V  e-(Y+ 5)r 

j=0  °  0  j=0  ° 

with 

A  o  =  32  (a  P  Y  5  )^'‘  /  ( a  +  p)^.  (43) 

The  computed  value  is  J  =  .455049,  the  same  as  the  formula  value. 
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VL  Conclusion 

All  expansion  methods  necessarily  have  their  limitations;  in  our  case  the 
number  of  harmonics  required  depend  upon  the  magnitude  of  the  product  of  each 
orbital  screening  constant  and  the  displacement  of  the  orbital  from  the  origin.  17 
However,  when  one  relates  this  basic  consideration  to  realistic  physical  problems, 
the  constraint  is  not  as  severe  as  might  be  supposed.  In  most  of  our  examples,  we 
achieved  excellent  results  with  use  of  only  13  harmonics.  We  note  that  exchange 
integrals  become  very  small  and  may  be  discarded  when  displacement  distances 
or  screening  constants  become  large;  also,  when  the  charge  overlap  between 
orbitals  is  essentially  zero  its  Coulomb  interaction  may  be  calculated  on  the  basis 
of  multipoles.  1^ 

The  efficiency  of  this  generalizable  method  for  the  evaluation  of  STO 
multicenter  molecular  integrals  may  be  significantly  improved  by  better 
programming,  the  use  of  look-up  tables,  and  vector  processing. 
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Figure  1.  Location  of  orbitals  X^(l),  X^{2),  and  X^i2)  for  the 

determination  of  the  three-center  exchange  integrals. 


Figure  2.  The  nine  regions  of  space  separately  calculated  for  the  evaluation  of 
the  three-center  exchange  integral. 
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Table  I.  The  values  of  thirteen  harmonics  for  the  three-center  exchange 
integral  over  the  top  half  of  the  quadrant,  and  the  bottom  half. 


f 


t\ 


fv  >o 
o  o  7 
^  'O  K 
fv  'O  ® 
•o  rv  r«) 
O  '«»v 

(ll  «4 

O  o 


in  ^  OD 

CO  o  n 

'O  Ov  ^ 

b1  ®  © 

ro  r  t  ^ 

4M  ©  Ok 

o 

«0  «H  ^ 

^  ©  -o 

r4  -*  o 

b't  «.«  o 

o  o  © 

©  o  © 

O  o  O 

O  o  © 

•  «  • 

©  o  o 

•  •  • 

^  'O 
^  li-jo 

ri  «-(bi 
M  -,0 
o  o  o 
o  o  o 

O  OO 
O  oO 


osjo-ioODm-r  rsrs.fio^*o 
O  "T  ^  ri  o  ri  ^3  til  n  >0  r)  -* 
'OrvCOo  ^  r>«  tnw*<  o  ciO 
'Oriffly)  Qr«)^  ooo  ooo 
^•H^f^Ooo  oOo  ooo 
C'Jf^OoOOO  OOo  OoO 
222®  2®®  ooc  OOO 
oooo  Ooo  ooo  Ooo 


f  I 


I 


I  i  i  I  t 


I  I 


CO*^ 


Oo>(hrtQf*i-^ 
o  ^  w  *0  rt  -0  ffl  r  i  ^  f-)  !L  o 
nrsf^o  ^  'O'O  ri**^o  o  oo 
'Ors^ro  ^-^o  oOo  ooo 
^  MOf^»^OoO  oOo  OoO 

~  ^!r2®2®®®®oooo 
OOOo  Ooo  ooo  Ooo 
OOOOOOO  oOo  Ooo 


I 


(  I  (  I  I 


I 


1  t 


r  Cf  2^^  O  o  o.  fs ^ 

f-iCDO^r  Ort-^  M3CKO  NofN 
—  fvtii.oNo'.'v 

f-i  t*i  ®  r  j  ^  (K  fv  f*|  o  o 
f^rx^in  •OQ.n  ..oo  ooo 
^  o  C'^  r\  r  t  o  o  o  o  o  ooo 
fv^CJoOoO  oooooo 
ooooooo  oooooo 


tft 


a  C4  ^  t'^s  IP  o-  ^  ^  ^  ooo 
<  f  I  o-  b*)  ®  r  j  ^  o  o  Coo 
mr^^o  rioo  ooo  oXo 
jw  '0«onooo©  oooooo 

o  ooo  O  oo  5  o5  5  oo 

22222®2  «>®o  Ooo 
OoOO  OoO  oOo  ooo 

I  I  I  I  I  I  I  III  III 

inr4  2-<r!-'' 

ooO-«  -c  ei  ^  ^^o  Ooo 
riesr^oo'n'*  oOo  ooo 

^C^^nOoO  oOo  OoO 
t^T'lOoOoO  oOoOoO 
2°2®2®'’  OOo  OoO 
222°2®®  ooo  ooo 
oooo  OoO  oOo  OoO 

'  I  '  I  >  I  I  III  I  I*  I 
22, -oTri  o-rj  — 


noM^-^ri'*  n«,o^v 
ncKOio  T-iin  s)Cin  nn  — 

222~'2'''’2  o”-<  ooo 
OQb^^  OOo  OoO 

^oN®r/^c  oooooo 
rviNtio  ooo  ooo  ooo 
OOOoOoO  oOo  Coo 


^bTffl  rs.rj^  o»^ri 
f^bTf-i  ©oo 
CO  ^  fnc^  rifirt  wO©  ©oo 
•^•COfv  ^  «4©  ©c©  ©oo 

O  •Ot>*®^OoOoCo©©© 
Cl^  f'lWOO©OOoO©©oO 
^O©©  ©©©  oC©  ©oO 
0©0©OoO  ©©o  ©oO 

*  I  *  »  I  I  <  I  I  I  I  I  I 

Ofs^  ©.nn 

oa^l*)<^  MS^tn  ©oo 

lotn^  o-  ©rj  o©©  ©oo 

©  «o  ('I  T'l  ri  o  o  o  ©  o  o  o  o 
'Oin'O^ooo  ©©o  ooo 
'em©©  ©o©  ©Oo  o  o® 
^oo©©©©  o©o©oO 
oo©o  o©©  oO©  ©oo 

I  '  I  <  I  I 


W 


I 


—  o  2®  5  o-  m  o  — tv  I'  rio- 
<  tr  » r  I  o  o,  r  /  o. ,'  _  n  ,0  f  * 
Ort»-nO‘ivl<^o.rvnnoO 
onoo  -.Oo  OoO 

M.  W— in-.0  oOo  OnO 

®  -JM  2  ®  ®  O  O  oo 
®l2o  Ooo  ooo  o5o 
nOOo  OoO  oOo  Oo© 


p 

o 

o 

fl 

IV 


V 

u 

V 

V 

u 


tu 

JS  . 

.-J  M 

c 

u  o 

^  ‘Sa 

V 

w  ^ 

•-? 

c  £ 

o  a 

£  to 

rt  "« 

■“  c 
c  3 

V  o 

V  o. 

•n  "a 

^  •a® 

'S  a 

I- 

<0  V 


V 


> 

o 

?  "5 

>  I- 
&0 
a;  o 

E^.S 


0) 

Iq 

CO 


0»«r#r»  ♦fO'O  N©».  0»<N 


